libcommonfunc.py

Go to the documentation of this file.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

## @file

##
# @package libcommonfunc Library of common functions for microwave filters synthesis
# @author J.M. Rius, J. Mateu, J.M. Tamayo, C. Collado, A. Padilla and J. O'Callaghan
# <p>Dpt. Signal Tehory and Communications, Universitat Politècnica de Catalunya (UPC).
# @version 2.0.3
# @date June 11, 2013
# <p>Acknowledgement: The software was developed in the frame of contract 21398/08/NL/GLC with the European Space Agency (ESA). Technical Offer was Christoph Ernst. 
# Further features were developed under contract UPC-C7767 with Thales Alenia Space España (TAS-E).
# Contributions to the definition of the software functionality and testing have been made by Christoph Ernst, Mónica Martínez Mendoza and other 
# ESA-ESTEC personnel, and Santiago Sobrino and Luis Roglá from TAS-E.
# <p>Copyright: &copy; Universitat Politècnica de Catalunya (UPC)  2009. <br>Contact: lossyfilters@tsc.upc.edu
# <p>This program is free software; you can redistribute it and/or modify it under the terms of
# the GNU General Public License as published by the Free Software Foundation; either 
# version 3 of the License, or (at your option) any later version.
# <p>This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; 
# without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
# <p>You should have received a copy of the GNU General Public License along with this program 
# in file "LICENSE.GPL3";  if not, download it from http://www.gnu.org/licenses/gpl-3.0.html  .
# <p>Additional permission under GNU GPL version 3 section 7:
# <p>If you modify this Program, or any covered work, by linking or combining it with the "Extra Filters Library" (libextrafilters), 
# the "Lossy Filters Library" (liblossyfilters) or the "License Check Library" (libchecklicense), or modified versions of that libraries, 
# containing parts covered by the terms described in files LICENSE.LIBEXTRAFILTERS and LICENSE.LIBLOSSYFILTERS, 
# the licensors of this Program grant you additional permission to convey the resulting work. 
# 

import re
import sys
import platform
import codecs
import os.path
import types
import csv
import numpy as np
import scipy as sc
import scipy.signal
import scipy.optimize as scop
from math import log10,  pi,  atan2,  sqrt,  floor
from PyQt4.QtGui import QMessageBox

from template import fileTemplate

# Disable 'Warning: divide by zero encountered in divide' in numpy. Result will be np. inf.
# These warning appears tipically when computing Q of resonators for a lossless filter.
np.seterr(divide='ignore')

#Multiprecision tests with mpmath
#import mpmath as mp
#mp.mp.dps = 100

#!!!!!!!!!!!!!!!!!!!!!!!!#
# Global variables:  #

##
# Dictionary of parameter keywords and variable names.
# This dictionary allows the developer to change a variable name in the code without having to change the corresponding keyword in all the existing parameter files.
parameterDict = { 'outName': 'outName',  'outDir': 'outDir',  'outDirName': 'outDirName',  'N': 'N',  'f0': 'f0',  'BW': 'BW',  
                            'filterTransferFunc': 'filterTransferFunc',  'chebyLr': 'chebyLr', 'qeZero': 'qeZero' ,  'qeLr':  'qeLr',  
                            'transImagZeros': 'transImagZeros',  'transComplexZeros': 'transComplexZeros', 'genChebyLr': 'genChebyLr',  
                            'genChebyLP': 'genChebyLP',  'LPalgor' : 'LPalgor', 'LPmaxIter' : 'LPmaxIter', 'LPNsamples' : 'LPNsamples', 
                            'LPcomplexZeros': 'LPcomplexZeros',  'LPrealZeros': 'LPrealZeros', 'LPfracBW': 'LPfracBW',  'LPripple': 'LPripple',  
                            'synthTech': 'synthTech',  'zerosArrang': 'zerosArrang',  'adapPredis': 'adapPredis',  'Qef': 'Qef',  'wFunc': 'wFunc',  
                            'finiteQo': 'finiteQo',  'Qo': 'Qo',  'QoPredis': 'QoPredis',  
                            'nuqK': 'nuqK',  'nuqK11c1': 'nuqK11c1',  'nuqK22c1': 'nuqK22c1',  'nuqK21c1': 'nuqK21c1',  'nuqK21c2': 'nuqK21c2',  'nuqK21c3': 'nuqK21c3',  
                            'nuqTech': 'nuqTech',  'nuqDelta': 'nuqDelta',  
                            'maxFreq': 'maxFreq',  'minFreq': 'minFreq',  'numFreq': 'numFreq' }

##    
# Number of significant digits to save in files
saveSignificantDigits = 14

##
# Number of significant digits to save in energy and power results.
# Usually less than saveSignificantDigits, because not much precision is needed.
saveSignificantDigitsEnergy = 4

##
# Module global variable equal either to None or to the GUI QMainWindow class instance, containing the application main window.
# It is equal to None (default) when the calling program is the console interface (mwfiltersgui.py --nogui).
# When creating the main window, the GUI executes the libcommonfunc.setGlobalVariablesLibCommonFunc() function in order to set this variable.
# The GUI methods do not use this variable. They access the mainWindow class through ParamEditDlg.mainWindow and MainWindow 'self' attributes.
mainWindow = None

##
# Module global variable equal to True if the GUI has been called with the -v flag command line option. Otherwise it is False.
verbose = False

##
# Module global variable equal to True if the GUI has been called with the --nogui flag command line option. Otherwise it is False.
nogui = False

##
# String containing console output generated before the GUI main window exists
preGuiOutput = ''

# End  global variables  #
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!#


##
# 
#     Parse error catch.
#     
class parseError(Exception): 
    pass    


##
# 
#     Synthesis error catch.
#     
class synthError(Exception): 
    pass    


##
# 
#     License error catch.
#     
class licenseError(Exception): 
    pass
    

#!!!!!!!!!!!!!!!!!!!!!!!!#
# Global functions:  #

_cfunc_complex = np.vectorize(complex)
    
#def mp2np(M):
#    """
#    Function to convert an mpmath complex matrix into a numpy array.
#    @param M = mpmath complex matrix
#    @returns npM = numpy array
#    """
#    npM = np.array( M.tolist() )
#    return _cfunc_complex(npM.ravel()).reshape(npM.shape)
#    
#    
#def np2mp(M):
#    """
#    Function to convert a numpy complex array or matrix into a mpmath matrix.
#    @param M = numpy complex array
#    @returns npM = mpmath matrix
#    """
#    return mp.matrix( M.tolist() )
    
    
##
# 
#     Function to allow the GUI set values to the global variables of libcomonfunc module.
#     This function is used in order to avoid circular imports.
#     @param variablesDict: Variable keyword argument list. Python creates a dictionary.
#     
def setGlobalVariablesLibCommonFunc(**variablesDict):    
    for key, value in variablesDict.items(): 
        globals()[key] = value
        
        
##
# 
#     Returns the value of the global verbose variable.
#     Useful to get the value from other modules that were loaded before the value of verbose was updated
#     
def readVerbose():
    return verbose
    

##
# 
#     Retrieves software package name and version from string.
#         @param st = String to parse for package name and version. It must be the docstring  (__doc__) of the calling module.
#         @return pkgname = Package name.
#         @return version = Package version.
#     
def pkgNameVersion(st):
    pkgname = re.search('@package (.*)', st).group(1)
    version = re.search('@version (.*)', st).group(1)
    return(pkgname, version) 
    

##
# 
#     Prints on standard output if console application.
#     or on log window if GUI application.
#         @param st = String.
#     
def myPrint(st):
    
    # Remove html tags from string:
    st = re.sub('<[^>]*?>','', st)
    
    if mainWindow is not None: 
        # In windows platform, remove '\r' end-of-line because Qt expects only \n
        if platform.system() in [ 'Windows',  'Microsoft',  'Vista' ]: st = st.replace('\r\n', '\n')
        
        # Another way to do the same:
        # if sys.platform.find('win') > -1: st = filter(lambda c: c!='\r', st )
        
        mainWindow.log_listWidget.addItem(unicode(st))     # addItem() method expects a string, not a Parameter class instance
        mainWindow.log_listWidget.setCurrentRow(mainWindow.log_listWidget.count()-1)    # Set log_listWidget current row to last one
        
    else: 
        # Check if console encoding is able to print non-ascii characters
        if hasattr(sys.stdout, 'encoding') and sys.stdout.encoding in [None, 'ascii']:      # If it is not
#            st = filter(lambda c: c<128, st )          # Remove non-ascii characters
            st = ''.join(map(lambda c: c if ord(c)<128 else '?', st ))     # Replace non-ascii characters by '?'
        print st


##
# 
#     Prints application header info.
#     The information shown is read from the file docstring.
#         @param applicationName = Name of application, usually os.path.basename(sys.argv[0]).
#         @param st = String to parse for package name and version. It must be the docstring  (__doc__) of the calling module.
#         @param importModule = Flag, set to True when this function is called to acknowledge successful module import.
#     
def printHeader(applicationName,  st,  importModule = False):
    
    if importModule:
        st = 'Successfully imported: %s; version %s\n' % pkgNameVersion(st) 
    else:
        st = u"\n%s\n%s" % (applicationName,  st)
        st = st.replace('@package',  'Package')
        st = st.replace('@author',  'Authors:')
        st = st.replace('@date',  'Date:')
        st = st.replace('@version',  'Version:')
        st = st.replace('<p>',  '\n')
        st = st.replace('<br>',  '\n')
        st = st.replace('&copy;',  '(c)')
        st = """%s \nPython %s on %s \n""" % (st,  platform.python_version(),  platform.system()) 
   
    myPrint(st)
    
    # Append st to  preGuiOutput if the GUI will be opened later
    if mainWindow is None and not nogui: 
        global preGuiOutput
        preGuiOutput = preGuiOutput + st


##
# 
#     Read parameters file and return a Parameter class instance.
#         @param fileName = Parameter file name.
#         @return Parameter class instance.
#     
def readParamFile(fileName):
    P = None 
    
    
    try:
        if fileName:
            for codec in ['utf-8',  'iso8859_15']:
                try:
                    myPrint('Openning file: %s' % os.path.basename(fileName))
                    file = codecs.open(fileName,  'r',  codec)
                except IOError,  err:
                    raise parseError,  'I/O error: %s' % err
                
                try:
                    stList = file.readlines()
                    break
                except UnicodeDecodeError, err:
                    #raise parseError, "Cannot open file with encoding '%s'<p>%s" % (codec,  err)
                    myPrint("Cannot open file with encoding '%s'" % (codec))
                
            myPrint("File sucessfully opened with encoding '%s'" % (codec)) 
                
            file.close()
    
            myPrint('Reading parameters from file: %s' % os.path.basename(fileName))
            P = Parameters(stList)

    except parseError,  err:
        criticalErrorMsg( 'ERROR in parameter file.<p>' + unicode(err) )
            
    return P


##
# 
#     Show critical error message. Prints error and exists in console mode. Displays QMessageBox in GUI mode.
#         @param st = Error string. It can contain <p> and <br> html tags, that are discarded in console mode.
#     
def criticalErrorMsg(st):
    if mainWindow is not None: 
        QMessageBox.critical(mainWindow,  "%s" % mainWindow.applicationName, st)
        
##@todo: See if MainWindow.cleanup() method can be used here.        
        # Close the [S] parameters plot and set all results to None
        if mainWindow.SPlotMagPhase is not None and not mainWindow.SPlotMagPhase.closed: 
            mainWindow.SPlotMagPhase.close()
            mainWindow.SPlotMagPhase = None
        if mainWindow.SPlotError is not None and not mainWindow.SPlotError.closed: 
            mainWindow.SPlotError.close()
            mainWindow.SPlotError = None
        CP = CM =  SP = None
    else:  
        st = st.replace('<p>',  '\n')
        st = st.replace('<br>',  ' ')
        st = re.sub('<[^>]*?>','', st)
        print st
        sys.exit(2)


##
# 
#     Show warning message. Prints warning in console mode. Displays QMessageBox in GUI mode.
#         @param st = Warning string. It can contain <p> and <br> html tags, that are discarded in console mode.
#     
def warningMsg(st):
    if mainWindow is not None: 
        QMessageBox.warning(mainWindow,  "%s" % mainWindow.applicationName, st)
    else: 
        st = st.replace('<p>',  '\n')
        st = st.replace('<br>',  ' ')
        st = st.replace('<b>',  '* ')
        st = st.replace('</b>',  ' *')
        print st


##
# 
#     Show information message. Prints message in console mode. Displays QMessageBox in GUI mode.
#         @param st = Information string. It can contain <p> and <br> html tags, that are discarded in console mode.
#     
def informationMsg(st):
    if mainWindow is not None: 
        QMessageBox.information(mainWindow,  "%s" % mainWindow.applicationName, st)
    else: 
        st = st.replace('<p>',  '\n')
        st = st.replace('<br>',  ' ')
        print st
        
        
##
# 
#     Show dialog to ask question.
#         @param st = Information string. It can contain <p> and <br> html tags, that are discarded in console mode.
#         @param buttons = QMessageBox buttons to show. (e.g.: QMessageBox.Yes | QMessageBox.No)
#         @return reply = button clicked ((e.g.: QMessageBox.Yes)
#     
def askQuestion(st,  buttons):
    if mainWindow is not None: 
        return QMessageBox.question(mainWindow,  "%s" % mainWindow.applicationName, st, buttons)
    else: 
        pass


##
# 
#     Nicely convert list of floats to string, rounding to 'prec' significant digits
#         @param listX = List of floats.
#         @param prec = Number of significant digits (int). Default 3.
#         @return st = String representation of x.
#     
def listStr(listX,  prec=3):
    st = '['
    for x in listX: st += '%.*g, ' % (prec,  x)
    return st[0:len(st)-2] + ']'


##
# 
#     Nicely convert complex number to string, rounding real and imaginary parts to 'prec' significant digits
#         @param x = Complex number.
#         @param prec = Number of significant digits (int).
#         @param format = None | 'DB' | 'MA' | 'RI' (string). Default is None. 
#         @param minWidth = Minimum field field (int). Default is 0.
#         @param eng = Use engineering notation in formats None and 'RI'. Default is False.
#                      Formatted strings of the form 0.x are allowed for exponents 0 and 9.
#         @param listFormat = Flag to return a list of strings instead of a string. The list contains [ Magnitude, Angle] or [Real, imag]. Useful to save in .csv format.
#                      
#         The format is borrowed from TouchStone s2p [S] parameter file definition.
#             None: Complex number, like 2.4-j5.3.
#             DB: Decibel Magnitude/Degree Angle.
#             MA: Linear Magnitude/Degree Angle (polar form).
#             RI: Real Part/Imaginary Part (rectangular form).
#         @return st = String representation of x.
#     
def complexStr(x,  prec,  format=None,  minWidth=0,  eng = False,  listFormat = False):

    re,  im = np.real(x),  np.imag(x)

    # Engineering notation
    if eng: 
        if re != np.inf:
            reExp = int(3*floor(log10(abs(re))/3)) if re != 0 else 0
            re = re / pow(10, reExp)
            if reExp in [-3,  6]: 
                if re >= 100:
                    re /= 1e3
                    reExp += 3
            stReExp = 'e' + str(reExp) if reExp != 0 else ''
        else: stReExp = ''
        
        if im != np.inf:
            imExp = int(3*floor(log10(abs(im))/3)) if im != 0 else 0
            im = im / pow(10, imExp)
            if imExp in [-3,  6]:
                if im >= 100:
                    im /= 1e3
                    imExp += 3
            stImExp = 'e' + str(imExp) if imExp != 0 else ''
        else: stImExp = ''
    else:
        stReExp, stImExp = '',  ''

    reim = abs(re/im) if abs(im) > 0 else np.inf
    imre = abs(im/re) if abs(re) > 0 else np.inf
        
    if format == None:
        st1 = '%+.*g%s' % (prec,  re,  stReExp)  if min(abs(re), reim) > 10**-prec else ''
        st2 = '%+.*g%sj' % (prec,  im,  stImExp) if min(abs(im), imre) > 10**-prec else ''
    elif format == 'RI':
        st1 = '%-+*.*g%s'   % (minWidth, prec,  re,  stReExp)  if min(abs(re), reim) > 10**-prec else '%-*s' % (minWidth,  '0')
        st2 = ' %-+*.*g%s' % (minWidth, prec,  im,  stImExp)  if min(abs(im), imre) > 10**-prec else ' %-*s' % (minWidth,  '0')
        st2 = st2.replace('+',  ' ')
    elif format == 'MA':
        mag,  angle = abs(x),  (180/pi)*atan2(np.imag(x),  np.real(x))
        st1 = '%-+*.*g'  % (minWidth, prec,  mag)     if abs(mag) > 10**-prec else '%-*s' % (minWidth,  '0')
        st2 = ' %-+*.*g' % (minWidth, prec,  angle)  if abs(angle) > 10**-prec else ' %-*s' % (minWidth,  '0')
    elif format == 'DB':
        if x==0:  # Set zero to -1000 dB
            mag = -1000
            angle = 0
        else: 
            mag,  angle = 20*log10(abs(x)),  (180/pi)*atan2(np.imag(x),  np.real(x))
        st1 = '%-+*.*g'  % (minWidth, prec,  mag)     if abs(mag)   > 10**-prec else '%-*s' % (minWidth,  '0')
        st2 = ' %-+*.*g' % (minWidth, prec,  angle)  if abs(angle) > 10**-prec else ' %-*s' % (minWidth,  '0')

    if listFormat == False:
        st = '%-*s' % (minWidth,  st1 + st2) 
        if len(st)==0: st = '0'
        if st[0]=='+': st = ' ' + st[1:len(st)]
    else:
        st = [st1,  st2]
    
    return st

# End global functions
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!#


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin class: Parameters

##
# 
#     Filter, synthesis and sweep parameters.
#     
class Parameters(object):
    
    ##
    # 
    #         Parse int, float, complex or string parameters from string stValue.
    #             @param key  = string containing parameter keyword.
    #             @param type = string containing parameter type: 'int', 'float', 'complex' or 'string'. 'Float' type allows a value of Inf, inf or INF.
    #             @param stValue = string containing parameter or a list of comma-separated parameters.
    #             @n Complex number format is 1+2j.
    #             @return parameter value.
    #         
    def readParam(self, key,  type,  stValue):

        # Remove trailing newline
        if stValue.endswith('\n'): stValue = stValue[:-1]
        
        # Check if parametr is empty
        if stValue.lstrip() == '': return None

        # Split the parameter values list
        stValue = stValue.split(',')
        
        # Parse stValue
        try:
            if    type == 'int':        n = [ int(elem) for elem in stValue ] 
            elif type == 'float':    n = [ np.inf  if elem.lower().lstrip().rstrip()=='inf' else float(elem) for elem in stValue ] 
            elif type == 'complex': n = [ complex(elem) for elem in stValue ] 
            elif type == 'string':   n = [ unicode(elem.lstrip().rstrip(),  'utf-8') for elem in stValue ] 
            else: raise parseError,  "Bad [%s] parameter type: '%s' " % (key, type)
        except ValueError,  err: 
            raise parseError,  "Bad [%s] parameter format: '%s' " % (key,  err)
            
        # Convert one-element list to element type
        if len(n) == 1: n = n[0]
        
        return n
        
        
    ##
    # 
    #         Initialize parameters class with list of strings.
    #         Each string is a line of the parameter file.
    #         For backwards compatibility of parameter files, missing fields are set to None.
    #             @param stList = List of strings. Each element of list is a line of input file.
    #                                         @n If stList is equal to None, an empty instance of Parameters class is created.
    #                                         (useful to create a new parameter file and for the copy funcion).
    #             @return Parameter class instance.
    #         
    def __init__(self, stList):
        if stList is not None:  # Create Parameter class instance from parameter information in stList
        
            # Extract (keyword, type, value)  from lines containing parameter values
            validLines = [ st for st in stList if st.find('=') >=0 and not st.startswith('#') ] 

            pat = re.compile('^.*\[(.+)\].*\{(.+)\}.*=(.+)$')
            matchList = [ pat.match(line) for line in validLines ]
            
            if matchList[0] is None: # Chebash file format
                myPrint("Reading Chebash file format")
                pat = re.compile('^(.+)=(.+)$')
                matchList = [ pat.match(line) for line in validLines ]
                chebashFormat = True
                
            else: # LossyFilters file format
                myPrint("Reading LossyFilters file format")
                chebashFormat = False
                
            # List of tuples containing (keyword, type, value) for each parameter
            paramList = [ match.groups(1) for match in matchList if match is not None]

            if chebashFormat:  # Read paramList in Chebash fornat
                self.readChebash(paramList)
            else:  # Read paramList in Lossyfilters fornat
                try:
                    for (key,  type, value) in paramList:
                        # Double backslashes have been removed by parameter matching, and must be restored. Otherwise, a utf8 decoding error will result from \nnn secuences
                        # It is necessary also to strip \r when reading windows files.
                        value = value.replace('\\', '\\\\').strip()     
                        command = "self.%s = self.readParam('%s', '%s', '%s')" % ( parameterDict[key], key,  type, value )
                        exec command
                except KeyError: raise parseError, u"Invalid keyword: '%s'" % key
            
            # For backwards compatibility of parameter files, set missing fields to None.
            for field in parameterDict.values():
                if field not in dir(self): exec "self.%s = None" % field
                
            #For backwards compatibility with Parameter files prior to implementation of Linear Phase optimization parameters
            if self.LPalgor is None:         self.LPalgor         = 'SLSQP'
            if self.LPmaxIter is None:      self.LPmaxIter     = 1000
            if self.LPNsamples is None:  self.LPNsamples  = 100
                
            # For backwards compatibility with Parameter files prior to implementation of asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11):
            if self.nuqK11c1 is None: self.nuqK11c1 = self.nuqK
            if self.nuqK22c1 is None: self.nuqK22c1 = self.nuqK
            if self.nuqK21c1 is None: self.nuqK21c1 = self.nuqK
            if self.nuqK21c2 is None: self.nuqK21c2 = self.nuqK
            if self.nuqK21c3 is None: self.nuqK21c3 = self.nuqK
            
            # Check if parameters have invalid values and make some basic parameter processing
            self.validate()
            
            ##
            # Flag that indicates if the class instance contains information.
            # empty==False or empty==True
            self.empty = False    
        
        else: # Create empty Parameter class instance with all fields set to None
            for key in parameterDict.values(): 
                command = "self.%s = None" %  key
                exec command
            self.empty = True
        # Set lists to empty list:
            self.transImagZeros = []
            self.transComplexZeros = []
            self.LPcomplexZeros = []
            self.LPrealZeros = []
            self.wFunc = []
            
            # Set values no to crash the GUI
            self.N = 2
            self.zerosArrang = 3
            
        ##
        # Path of output files, including directory and file name
        self.outDirName = None
        
    
    ##
    # 
    #         Translate a list of Chebash parameters into a list of equivalent lossyfilters parameters
    #         @param paramList = List of tuples (Chebash_Parameter_Name, value) where value is a string
    #         @return paramList = List of tuples (LossyFilters_Parameter_Name, type, value) where type and value are strings
    #         
    def readChebash(self, paramList):
        
##        def normRealFreq(f):
##            """
##            Auxiliary function
##            """
##            return (f/self.f0 - self.f0/f)*self.f0/self.BW

        
##        for (keyChebash,  stValue) in paramList:

        # We have to assume that Chebash parameters  appear in the same line order.
        # We cannot use keywords to find the parameters, since there are repeated keywords (e.g. 'real', 'imag')
        if paramList[0][0] == 'n': self.N = int(paramList[0][1]) 
        else: raise parseError,  u"parameter 'n' missing or misplaced in Chebash file"

        if paramList[1][0] == 'fo': self.f0 = float(paramList[1][1])  * 1e6
        else: raise parseError,  u"parameter 'fo' missing or misplaced in Chebash file"

        if paramList[2][0] == 'Qo': self.Qo = float(paramList[2][1])
        else: raise parseError,  u"parameter 'Qo' missing or misplaced in Chebash file"

        if paramList[3][0] == 'wp': self.BW = float(paramList[3][1])  * 1e6
        else: raise parseError,  u"parameter 'wp' missing or misplaced in Chebash file"

        if paramList[4][0] == 'ro': self.genChebyLr = float(paramList[4][1])
        else: raise parseError,  u"parameter 'ro' missing or misplaced in Chebash file"

        if paramList[5][0] == 'puntos': self.numFreq = int(paramList[5][1]) 
        else: raise parseError,  u"parameter 'puntos' missing or misplaced in Chebash file"

        if paramList[6][0] == 'n_ceros': n_ceros = int(paramList[6][1]) 
        else: raise parseError,  u"parameter 'n_ceros' missing or misplaced in Chebash file"
        
        if n_ceros > 0:
            if paramList[7][0] == 'ceros': self.transImagZeros = [ self.f0 + float(st)*1e6 for st in paramList[7][1].split(',') ]
            else: raise parseError,  u"parameter 'ceros' missing or misplaced in Chebash file"       
            if len(self.transImagZeros) != n_ceros: raise parseError,  u"Incorrect number of transmission zeros in Chebash file"
        else:
            self.transImagZeros = []
            
        if paramList[8][0] == 'n_ceros_eq': n_ceros_eq = int(paramList[8][1]) 
        else: raise parseError,  u"parameter 'n_ceros_eq' missing or misplaced in Chebash file"        

        if n_ceros_eq > 0:
##            if paramList[9][0] == 'real': realEqZs   = [ normRealFreq(self.f0 + float(st)*1e6) for st in paramList[9][1].split(',') ]
            if paramList[9][0] == 'real': realEqZs   = [ float(st) for st in paramList[9][1].split(',') ]
            else: raise parseError,  u"parameter 'real' (for equalization ceros) missing or misplaced in Chebash file"       
            
            if paramList[10][0] == 'imag': imagEqZs   = [ self.f0 + float(st)*1e6 for st in paramList[10][1].split(',') ]
            else: raise parseError,  u"parameter 'imag' (for equalization ceros) missing or misplaced in Chebash file"       
            
            self.transComplexZeros = [ complex(val[0],val[1]) for val in zip(realEqZs, imagEqZs) ]
            if len(self.transComplexZeros) != n_ceros_eq: raise parseError,  u"Incorrect number of equalization zeros in Chebash file"
        else:
            self.transComplexZeros = []
            
        # Remove complex transmission zeros with symmetric real part
        if n_ceros_eq > 0:
            nonSymmZeros = [ ]      # List of zeros removing symmetric in the real part
            for TZ in self.transComplexZeros:
                if -TZ.conjugate() not in nonSymmZeros: nonSymmZeros.append(TZ)
            self.transComplexZeros = nonSymmZeros            
            
        # Default values for variables not in Chebash files
        self.outName = 'Chebash'
        self.outDir = '.'
        self.filterTransferFunc = 'Generalized Chebyshev'
        self.genChebyLP = 0
        self.synthTech = 'Minimum Insertion Loss'
        self.finiteQo = 1
        self.nuqTech = 'kS21+kS11' 
        self.minFreq = self.f0 - 3*self.BW
        self.maxFreq = self.f0 + 3*self.BW
        
##        ##
##        # List of Chebash parameters to ignore
##        chebashIgnoreList = ['n_ceros', 'n_frecpropias_eq', 'real', 'imag', 'tipo_transformacion', 'graf[1]', 'graf[2]', 'graf[3]', 'graf[5]', 'graf[6]', 'graf[7]', 
##                                            'graf[8]', 'graf[9]', 'graf[10]', 'graf[11]', 'graf[12]', 'graf[13]', 'graf[14]', 'graf[15]', 'graf[16]', 'graf[17]', 'graf[18]']
##                                            
##        ##
##        # List of Chebash Mask parameters not to read
##        chebashMaskList = [ 'n_at', 'f_at', 'A_at', 't_at', 'offsetmba', 't_offsetmba', 'n_perd', 'f_abp', 'a_abp', 't_abp', 'offsetmbp', 't_offsetmbp', 
##        'nr_retardo', 'f_retardo', 'r_retardo', 't_retardo', 'offsetrg', 't_offsetrg', 'nptos_pendht', 'f_pendht', 'd_pendht', 't_pendht', 'offsetpendht', 't_offsetpendht',
##        'nptos_pendre', 'f_pendre', 'd_pendre', 't_pendre', 'offsetpendre', 't_offsetpendre', 'nptos_gdst', 'f_gdst', 'v_gdst', 't_gdst', 'offsetgdst', 't_offsetgdst', 
##        'nptos_gdsp', 'f_gdsp', 'v_gdsp', 't_gdsp', 'offsetgdsp', 't_offsetgdsp', 'n_frecpropias_eq', 'variacion_temp', 'margen_temp', 'modificar_mascaras', 'margen_gdsp']

##        processingComplex = False   # Will be True when we start processing a 3-line complex parameter
##        for (keyChebash,  stValue) in paramList:
##            
##            if not processingComplex:    # int, float, string or 1st line of a complex parameter
##                try:
##                    (keyLossyFilters, type, scale) = chebashTranslateDict[keyChebash]
##                except KeyError, err: 
##                    if keyChebash not in chebashIgnoreList + chebashMaskList: myPrint('Unknown parameter key in Chebash file: %s' % err)
##                    continue
##                    
##                if type in ['int', 'float']:
##                                        
##                    # First process the special cases
##                    if keyChebash == 'graf[4]':
##                        fmin, fmax = tuple( stValue.split(',')[0:2] )
##                        paramListLF.append( ('minFreq', 'float', str(float(fmin)*scale) ) )
##                        paramListLF.append( ('maxFreq', 'float', str(float(fmax)*scale) ) )
##                        continue
##                    else:
##                        if stValue != 'X':  # In Chebash, 'X' indicates an empty parameter
##                            exec("value = [ scale* %s(n) for n in stValue.split(',') ]" % type)
##                            if len(value) == 1:     # Not a list
##                                value = str(value[0])
##                            else:   # List
##                                value = str(value).replace('[', '').replace(']', '')
##                        else:   # Empty parameter
##                            value = ' '
##                        
##                elif type == 'string':
##                        value = stValue
##                elif type == 'complex':
##                    processingComplex = True
##                    NComplex = int(stValue)
##                else:
##                    raise parseError, 'Invalid data type'
##                    
##            else:   # 2nd and 3rd lines of a complex parameter
##                if keyChebash == 'real':
##                    if NComplex > 0: 
##                        complexRe = [ float(n)*scale for n in stValue.split(',') ]
##                elif keyChebash == 'imag': 
##                    if NComplex > 0:
##                        complexIm = [ float(n)*scale for n in stValue.split(',') ]
##                        value = str([ complex(nReal, nImag) for nReal, nImag in zip(complexRe, complexIm) ]).replace('[', '').replace(']', '').replace('(', '').replace(')', '')
##                    else:    
##                        value = ' '
##                    processingComplex = False
##                    
##            # Add lossyfilters parameter, for normal or for special case
##            paramListLF.append((keyLossyFilters, type, value))
##            
##        return paramListLF + chebashDefaultList
        
        
    ##
    # 
    #         Function to check if parameters have invalid values and make some basic parameter processing.
    #         
    def validate(self):
        
        if not isinstance(self.outName, types.StringTypes) or len(self.outName) == 0:
            raise parseError,  u"Name of results output files is not a valid string"
            
# DO NOT CHECK OUTDIR KNOW. CKECK IT AFTER COMPUTATION, JUST BEFORE SAVING
# THIS WILL ALLOW MODIFICATION OF OUTDIR USING THE GUI
#
#        if self.outDir is not None: 
#            # Normalize and expand output directory name, just in case
#            self.outDir = os.path.normpath(self.outDir)
#            self.outDir = os.path.expanduser(self.outDir)
#            if not isinstance(self.outDir, types.StringTypes) or not os.path.isdir(self.outDir):
#                raise parseError,  u"Output directory parameter '%s' is not an existing directory name" % self.outDir
#                
#        # Create output name using outDir and outName
#        if self.outDir is not None: self.outDirName =  os.path.join(self.outDir,  self.outName)
#        else:                                   self.outDirName =  self.outName
#            
#        # Normalize output path just in case 
#        self.outDirName = os.path.normpath(self.outDirName)
        
        if self.filterTransferFunc.title() not in [ u'Butterworth', u'Chebyshev', u'Quasieliptic',  u'Generalized Chebyshev' ]:
            raise parseError,  u"Invalid filter topology type '%s'" % self.filterTransferFunc
        
        if self.synthTech.title()  not in [ u'Minimum Insertion Loss', u'Predistortion',  u'Prescribed Insertion Loss' ]:
            raise parseError,  u"Invalid synthesis technique '%s'" % self.synthTech
            
        if self.nuqTech not in [ u'kS21+kS11',  u'kS21',  u'kS21+pole/zero' ]:
            raise parseError,  u"Invalid Prescribed Insertion Loss technique = '%s'" % self.nuqTech
            
        if self.numFreq <= 0:
            raise parseError,  u"Sweep parameters: Number of samples is not > 0"
            
        if self.minFreq >= self.maxFreq:
            raise parseError,  u"Sweep parameters: Minimum frequency is not smaller than Maximum frequency"
            
        if self.minFreq <= 0:
            raise parseError,  u"Minimum frequency must be larger than 0 Hz"
            
        if self.f0<= 0:
            raise parseError,  u"Central frequency must be larger than 0"
            
        if self.BW<= 0:
            raise parseError,  u"Bandwidth must be larger than 0"
            
        if self.numFreq <= 1:
            raise parseError,  u"Number of frequency samples must be larger than 1"
        
        if self.filterTransferFunc.title() ==  u'Chebyshev' and self.chebyLr <= 0:    
            raise parseError,  u"Chebyshev filter return losses must be larger than 0"
            
        if self.filterTransferFunc.title() ==  u'Quasieliptic':
            if self.qeLr <= 0:      raise parseError,  u"Quasieliptic filter return losses must be larger than 0"   
            if self.qeZero <= 1: raise parseError,  u"Quasieliptic filter transmission zero must be larger than 1"

        # Lists that could be empty in the parameter file
        if self.transImagZeros is None: self.transImagZeros = [ ]
        if self.transComplexZeros is None: self.transComplexZeros = [ ]
        if self.LPcomplexZeros is None: self.LPcomplexZeros = [ ]
        if self.LPrealZeros is None: self.LPrealZeros = [ ]
        if self.wFunc is None: self.wFunc = [ ]
        
        # Be sure that lists of zeros are really lists
        if type(self.transImagZeros) != types.ListType:  self.transImagZeros =  [ self.transImagZeros ]
        if type(self.transComplexZeros)  != types.ListType:  self.transComplexZeros =   [ self.transComplexZeros ]
        if type(self.LPcomplexZeros) != types.ListType:  self.LPcomplexZeros =  [ self.LPcomplexZeros ]
        if type(self.LPrealZeros) != types.ListType:  self.LPrealZeros =  [ self.LPrealZeros ]
        
        # Lists that should be a one-element list when it is a number
        if type(self.LPcomplexZeros).__name__ in ['int',  'float',  'complex']: self.LPcomplexZeros = [ self.LPcomplexZeros ]
        if type(self.LPrealZeros).__name__        in ['int',  'float',  'complex']: self.LPrealZeros = [ self.LPrealZeros ]
 
        if self.filterTransferFunc.title() ==  u'Generalized Chebyshev':
            if self.genChebyLr <= 0:
                raise parseError,  u"Generalized Chebyshev filter return losses must be larger than 0"
            
            if self.genChebyLP and any( np.real(self.LPcomplexZeros + self.LPrealZeros)==0 ):
                raise parseError, 'Linear Phase optimization zeros cannot be purely imaginary'
                
            if self.genChebyLP and any( np.imag(self.LPcomplexZeros)<=0 ):
                raise parseError, 'Imaginary part of Linear Phase optimization complex zeros cannot be negative or 0 Hz'
                
            if len(self.transImagZeros)>0 and any(np.array(self.transImagZeros)<=0):
                raise parseError,  'Imaginary transmission zeros cannot be negative or 0 Hz'
                
            if len(self.transComplexZeros)>0 and any(np.imag(self.transComplexZeros) <=0):
                raise parseError,  'Imaginary part of complex transmission zeros cannot be negative or 0 Hz'

        if self.synthTech.title() == u'Minimum Insertion Loss':
            if self.finiteQo and self.Qo <=0:     raise parseError,  u"Q of implementation must be larger than 0"
            
        if self.synthTech.title() == u'Predistortion':
            if self.QoPredis <=0:     raise parseError,  u"Predistortion Q of implementation must be larger than 0"
            if self.Qef <=0:     raise parseError,  u"Predistortion Q to emulate must be larger than 0"
            
            if type(self.zerosArrang).__name__ == 'int' and self.zerosArrang>=1 and self.zerosArrang<=4:  pass
            elif type(self.zerosArrang).__name__ == 'list':
                if len(self.zerosArrang) != self.N: raise parseError, u"Length of Predistortion factorization method boolean list is different than filter order"
                if type(self.zerosArrang[0]).__name__ == 'bool': pass
                elif type(self.zerosArrang[0]).__name__ == 'int': 
                    for n in range(self.N): self.zerosArrang[n] = bool(self.zerosArrang[n])
                else: raise parseError, u"Predistortion factorization method is not a list of Booleans or integers than can be converted to Booleans"
                
            else: raise parseError, u"Predistortion factorization method is not an integer in the range [1,4] or a list of Booleans"
            
            if self.adapPredis == 1:
                if len(self.wFunc) != self.N: raise parseError,  u"Length of adaptive predistortion weights list different than filter order"
                
        if self.synthTech.title() == u'Prescribed Insertion Loss':
            if self.nuqTech == 'kS21+kS11':
                if self.nuqK11c1 <=0:     raise parseError,  u"Prescribed Insertion Loss k21+kS11: k11 must be larger than 0"
                if self.nuqK22c1 <=0:     raise parseError,  u"Prescribed Insertion Loss k21+kS11: k22 must be larger than 0"
                if self.nuqK21c1 <=0:     raise parseError,  u"Prescribed Insertion Loss k21+kS11: k21 must be larger than 0"
            if self.nuqTech == 'kS21':                
                if self.nuqK21c2 <=0:     raise parseError,  u"Prescribed Insertion Loss k21: k21 must be larger than 0"
            if self.nuqTech == 'kS21+pole/zero':                
                if self.nuqK21c3 <=0:     raise parseError,  u"Prescribed Insertion Loss kS21+pole/zero: k21 must be larger than 0"
                if self.nuqDelta <= 0: raise parseError,  u"Prescribed Insertion Loss kS21+pole/zero: Delta must be larger than 0"
                eps = 0.05 # Tolerance in the validation of Delta parameter
                minDelta = 20 * 10**(-self.nuqK21c3/10) - eps
                if self.nuqDelta < minDelta: warningMsg("WARNING: Prescribed Insertion Loss kS21+pole/zero:<p>Delta parameter smaller than the<br>minimum recommended value (%.2g)" % minDelta)


    ##
    # 
    #         Function to copy instances of the Parameters class (copies the self instance).
    #             @return Pcopy = New instance of Parameters class copied from the self instance.
    #         
    def copy(self):
        Pcopy = Parameters(None)
        
        if self.empty == False:  # If the class instance to copy is not empty, copy it. Otherwise just return the empty Pcopy
            for key in parameterDict.values(): 
                command = "Pcopy.%s = self.%s" %  (key,  key)
                exec command
            Pcopy.empty = False

        return Pcopy
    
    
    ##
    # 
    #         Special method to print parameter class (prints the self instance).
    #         Floating point variables are printed as strings using str() because str() is inmune to roundoff errors in the floating point representation.
    #         
    #             @return st = Human redable printing of the parameters (string).
    #         
    def __str__(self):
        if self.empty == True:
            print ''
            return
            
        # Pattern to find parameter keywords
        pat = re.compile('^.*\[(.+)\].*')
        
        out = ''
        for line in fileTemplate.split('\n'):
            # Substitue parameter values in template
            if line.find('=') >=0 and not line.startswith('#'):
                key = pat.match(line).groups(0)[0]
                value = eval('self.' + parameterDict[key])
                
                # Round floats and lists of floats (converted to string in order to represent properly the rounded number)
                if type(value) == types.FloatType or type(value) == types.ComplexType:
                    value = complexStr(value,  saveSignificantDigits, eng = True)
                elif type(value) == types.ListType:
                    value = ', '.join( [ complexStr(val, saveSignificantDigits,  eng = True) for val in value ] )
                
                if value is not None:
                    st = 'Inf' if value==np.inf else unicode(value).strip()
                else:                           
                    st = ''
                    
                if type(value) == types.UnicodeType: 
                    line = line + filter(lambda c: c not in "[]", st )         # Allow () in strings
                else:
                    line = line + filter(lambda c: c not in "()[]", st )      # Remove () from complex numbers
                                
            # Write platform-dependant end-of-line
            if platform.system() in [ 'Windows',  'Microsoft',  'Vista' ]:
                out = out + line + '\r\n'
            elif platform.system() in [ 'Linux' ]:
                out = out + line + '\n'
            else: # Mac platform
                out = out + line + '\r'
            
        return out
        
# End class:  Parameters   #
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin class: Sparameters

##
# 
#     [S] parameters class.    
#     
class Sparameters(object):
    
    ##
    # 
    #         Initialize frequency axis attributes and compute S parameters from Characteristic Polynomials.
    # 
    #         @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #         @param CP = CharPolys  class instance, containing characteristic polynomials and constants.
    #         @param FT = Frequency transform class instance, containing methods to normalize frequency and unnormalize group delay.
    #         @param freq = Frequency axis (Hz). If equal to None, the frequency sweep based on P.minFreq, P.maxFreq and P.numFreq will be used. Default None.
    #         @param checkPassive = If checkPassive is True, a warning will be raised when there are [S] parameters larger than 0 dB. Default True.
    #         
    #         [S] parameters computed from characteristic polynomials.
    #         \f[ 
    #             S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)} 
    #         \f]
    #         \f[
    #             S_{21} = S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)}
    #         \f]
    #         \f[
    #             S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)}
    #         \f]
    #         except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where it is computed as
    #         \f[
    #             S_{22} = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right)  + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11}
    #         \f]
    #         In this "Prescribed Insertion Loss technique" case 1, the correct lossy values of \f$  S_{11} \f$ and \f$  S_{21} \f$
    #         are obtained because  \f$ \epsilon_R \f$ and \f$ \epsilon \f$  have been divided by \f$ k_{11} \f$ and \f$ k_{21} \f$ , respectively,
    #         when processing the polinomials at liblossyfilters.nonUniformQ() function .
    #         
    def __init__(self, P,  CP,  FT, freq = None, checkPassive=True):
        if CP == None: return
        
        ##
        # Fractional bandwidth of the FT used to compute this [S] parameters.
        self.FBW = FT.FBW
        
        ##
        # Central frequency of the FT used to compute this [S] parameters.
        self.f0 = FT.f0
        
        ##
        # Frequency transform used to compute this [S] parameters.
        self.FT = FT
        
        ##
        # Frequency axis (Hz)
        self.freq = None
        
        if freq is None: 
            self.freq = np.linspace(P.minFreq,  P.maxFreq,  P.numFreq)
        else:
            self.freq = freq
            
        ##
        # Normalized \f$ \omega ' \f$ axis
        self.omega = FT.normFreq(self.freq)
        
        # Value of quality factor used to compute dissipation factor sigma.
        # It is infnite in Prescribed Insertion Loss and in Minimum Insertion Loss ideal filter, otherwise it is the value in parameters file.
        if P.synthTech == u'Minimum Insertion Loss':
            if P.finiteQo: Qo = P.Qo
            else: Qo = np.Inf
        elif P.synthTech == u'Predistortion': 
            Qo = P.QoPredis
        elif P.synthTech == u'Prescribed Insertion Loss': 
            Qo = np.Inf
        
        ##
        # Dissipation factor. [Cameron eq. (3.118)] [TN  101.1 eq. (14)].
        #
        # \f$ \sigma = \frac{f_0}{\Delta f}  \frac{1}{Q_0} \f$
        #
        self.sigma = 1 / (Qo * FT.FBW)
        
        ##
        # Complex frequency variable in normalized s-plane.
        # \f$ s = \sigma + j \omega ' \f$ 
        self.s = self.sigma + 1j * self.omega
        
        ##
        # Parameter \f$ S_{11} \f$ computed from characteristic polynomials.
        # \f$ S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)} \f$
        self.S11 = None
        
        ##
        # Parameter \f$ S_{12} \f$ computed from characteristic polynomials.
        # \f$ S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)} \f$
        self.S12 = None
        
        ##
        # Parameter \f$ S_{21} \f$ computed from characteristic polynomials.
        # \f$ S_{21} = \frac{1}{\epsilon} \frac{P(s)}{E(s)} \f$
        self.S21 = None
        
        ##
        # Parameter \f$ S_{22} \f$ computed from characteristic polynomials.
        # \f$ S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)} \f$ 
        # [Cameron eq. (6.15) pp. 212]
        #
        # except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where it is computed as
        # \f$ S_{22} =  1 - \left ( \frac{k_{21}}{k_{11}} \right) ^2 + \frac{k^2_{21}}{k_{11}k_{22}} S_{11_lossless} \f$
        #
        self.S22 = None
        
        # Compute self.S11,self.S12, self.S21, self.S22
        self.fromCharPolys(P, self.s, CP.Ps, CP.Es, CP.Fs, CP.eps, CP.epsR)
        
        if checkPassive:
            epstol = 1e-3
            S11max,  S21max,  S22max = max(abs(self.S11)), max(abs(self.S21)), max(abs(self.S22))
            if max([S11max,  S21max,  S22max]) > 1+epstol: 
                warningMsg("""WARNING: [S] parameters computed from CP are larger than 0 dB at some frequency points<p>
                        max[S11] = %.2f dB, max[S21] = %.2f dB, max[S22] = %.2f dB
                      """ % (20*log10(S11max), 20*log10(S21max), 20*log10(S22max)) )
#            if any(abs(self.S11) > 1+epstol) or any(abs(self.S21) > 1+epstol) or any(abs(self.S22) > 1+epstol): 
#                warningMsg("WARNING: [S] parameters computed from CP are larger than 0 dB at some frequency points")
        
        ##
        # Group delay for the LowPass prototype
        self.groupDelayLP = self.groupDelayfunc(CP.Es,  CP.Ps,  self.s)
        # Computation of group delay as a derivative of the phase of S21
        # Fails when sampling S21 near the transmission zeros
        #self.groupDelayLP = - np.ediff1d(self.phaseS21) / np.ediff1d(self.omega)
        # Normalized \f$ \omega \f$ at the sampling points of the group delay
        #self.omegaGD = 0.5*(self.omega[1:] + self.omega[:-1])
        
        ##
        # Group delay for the unnormalized filter (sec)
        self.groupDelay = FT.unormGroupDelay(self.freq,  self.groupDelayLP)
        self.groupDelay[ self.groupDelay<0 ] = 0 # CE prefers that the GD goes negative but TAS-E prefers truncation of negative values
        
        ##
        # Phase of transfer function (rad)
        self.phaseS21 = np.unwrap(np.angle(self.S21))
        
        ##
        # Phase of return losses (rad)
        self.phaseS11 = np.unwrap(np.angle(self.S11))
        
        ##
        # Phase of inverse return losses (rad)
        self.phaseS22 = np.unwrap(np.angle(self.S22))
        
        # Best fit to linear phase straight line, fitting to straight line in the whole band
        #coef = np.polyfit(self.omega,  self.phaseS21,  1)   
        
        if freq is None:
            # Best fit to linear phase straight line, fitting to straight line in the center of the band in an interval equal to half the bandwidth
            NBW =  int( P.numFreq * P.BW/(P.maxFreq - P.minFreq) )  # Number of samples corresponding to BW
            N1 = int(P.numFreq/2 - NBW/4)-1
            N2 = int(P.numFreq/2 + NBW/4)+1
            coef = np.polyfit(self.omega[N1:N2],  self.phaseS21[N1:N2],  1)  
        
            ##
            # Deviation from linear phase.
            # It is the error between the phase of S21 and the straight line that best fits the phase in the center of the band in an interval equal to half the bandwidth.
            self.deviationLP = self.phaseS21 - (self.omega*coef[0] + coef[1])
        
        ##
        # Parameter \f$ S_{11} \f$ computed from a coupling matrix.
        self.S11M = None
        
        ##
        # Parameter \f$ S_{21} \f$ computed from a coupling matrix.
        self.S21M = None
        
        ##
        # Parameter \f$ S_{22} \f$ computed from a coupling matrix.
        self.S22M = None
        
        ##
        # Power dissipated at resonators.
        # 2-D numpy array of shape (Nf,  Nr), where Nf is the number of frequencies and Nr the number of resonators.
        self.powerRes = None
        
        ##
        # Energy stored at resonators.
        # 2-D numpy array of shape (Nf,  Nr), where Nf is the number of frequencies and Nr the number of resonators.
        self.energyRes = None
        
        ##
        # Power dissipated at resistive couplings.
        # Tuple ( (m, n), Gmn, powerDirect, powerReverse ) where power is a 1-D numpy array with one element for each frequency.
        # Only entries n>m are stored.
        self.powerCoup = None
        
        ##
        # Energy stored at magnetic or electric couplings.
        # Tuple ( (m, n), Bmn, energyDirect, energyReverse ) where energy is a 1-D numpy array with one element for each frequency.
        # Only entries n>m are stored.
        self.energyCoup = None
        
        ##
        # Power dissipated at resonators, with reverse excitation.
        # 2-D numpy array of shape (Nf,  Nr), where Nf is the number of frequencies and Nr the number of resonators.
        self.powerResRev = None
        
        ##
        # Energy stored at resonators, with reverse excitation.
        # 2-D numpy array of shape (Nf,  Nr), where Nf is the number of frequencies and Nr the number of resonators.
        self.energyResRev = None
    
    
    ##
    # 
    #         Compute S-Parameters (\f$S_{11}\f$, \f$S_{12}\f$, \f$S_{21}\f$, and \f$S_{22}\f$) from the characteristic polynomials.
    #         Evaluation points are taken from attribute Sparameters.s . 
    #         This method will be called by the BW and TZs optimization routines.
    #         Output in atributes Sparameters.S11, Sparameters.S12, Sparameters.S21 and Sparameters.S22 .
    #         
    #         @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #         @param evalS = Vector with the normalized complex frequencies (s) where we want to compute the group delay.
    #         @param Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$.
    #         @param Es = Characteristic polynomial of the denominator of \f$S_{11}\f$ and \f$S_{21}\f$.
    #         @param Fs = Characteristic polynomial of the numerator of \f$S_{11}\f$.
    #         @param eps = Characteristic constant \f$ \epsilon \f$ in the denominator of \f$ S_{12} \f$ and \f$ S_{21} \f$
    #         @param epsR = Characteristic constant \f$ \epsilon_R \f$ in the denominator of \f$ S_{11} \f$
    #         
    def fromCharPolys(self, P, evalS, Ps, Es, Fs, eps, epsR):
        
        # Evaluate polynomials for s: The result is of type numpy.array
        aFs = Fs(evalS)
        aPs = Ps(evalS)
        aEs = Es(evalS)
        
        self.S11 = (1/epsR) * aFs / aEs        
        self.S12 = (1/eps)   * aPs / aEs        
        self.S21 = self.S12        
        
        if P.synthTech.title() == u'Prescribed Insertion Loss' and P.nuqTech == u'kS21+kS11':
            k11 = 10**(-P.nuqK11c1/20.0)  # Insertion lossess
            k22 = 10**(-P.nuqK22c1/20.0)  # Insertion lossess
            k21 = 10**(-P.nuqK21c1/20.0)  # Insertion lossess
            self.S22 = k22 * ( 1 - (k21/k11)**2 ) + ( k21/k11)**2 * (self.S11)
        else:
            self.S22 = (1/epsR) * (-1)**P.N * aFs.conj() / aEs
            # Do not use F22s: S22 =(1/epsR) * F22s / Es , because libfreefilters.polyConj() is ony for 's' in the imaginary axis.    
    
    
    ##
    # 
    #         This function returns the group delay as defined in [Cameron eq. 13.7,  pp. 475-476] evaluated in the normalized complex frequency vector evalS .
    #         @param Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
    #         @param Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$.
    #         @param evalS = Vector with the normalized complex frequencies (s) where we want to compute the group delay.
    #         @return groupD = Vector with the group delay at each frequency.
    #         The group delay is computed from the formula:
    #         \f[
    #         \tau_{i}= Re\left[ \frac{E'(s_{i})}{E(s_{i})} -  \frac{P'(s_{i})}{P(s_{i})} \right]
    #         \f]
    #         where \f$s_{i}\f$ are the frequency samples.
    #         
    def groupDelayfunc(self, Es,  Ps, evalS):
        groupD = np.real( Es.deriv()(evalS)/Es(evalS) - Ps.deriv()(evalS)/Ps(evalS))
        return groupD        


    ##
    # 
    #         Compute S-Parameters (\f$S_{21}\f$ and \f$S_{11}\f$) from a coupling matrix.
    #         Evaluation points are taken from attribute Sparameters.s .
    #         Output in atributes Sparameters.S11M and Sparameters.S21M .
    #         @param MatQ = MatrixQ class instance, containing the coupling matrix to process.
    #         @param evalS = Normalized frequency axis (imaginary). If equal to None, it is set to 1j * np.imag(self.s). Default None.
    #         @param checkPassive = If checkPassive is True, a warning will be raised when there are [S] parameters larger than 0 dB. Default True.
    #         @return stErrorCM = String with a message about the error in [S] computed from M coupling matrix versus self.S11 and self.S21.
    #         @return CM_error = Relative error in [S] computed from M coupling matrix versus self.S11 and self.S21 (float).
    #         
    #         By definition, \f$S_{21}(s)\f$ and \f$S_{11}(s)\f$ can be found as:
    #         \f[
    #             S_{21}(s) = 2\sqrt{R_{S}R_{L}} \, Z'_{N+2,1}(s)
    #         \f]
    #         \f[
    #             S_{11}(s) = 1 - 2R_{S}Z'_{1,1}(s)
    #         \f]
    #         \f[
    #             S_{22}(s) = 1 - 2R_{L}Z'_{N+2,N+2}(s)
    #         \f]
    #         where:
    #         \f[ Z'(s)=Y'(s)^{-1} \f]
    #         \f[ Y'(s)=jM+sI+G \f]
    #         with \f$I \f$ a diagonal matrix with the node scaling factor squared in the diagonal except in the non-resonant nodes positions, where it is zero. 
    #         It means putting zeros in the extraNodesS first positions of the diagonal and in the last extraNodesL positions.
    #         
    #         \f[ G = \begin{bmatrix} G_{S} & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & & G_{L}\end{bmatrix} \f]
    #         @todo Find out why sign of S11 from CM is opposite than from CP. The sign is changed below to force agreement with S11 from CP.
    #         
    def fromCouplingMatrix(self,  MatQ, evalS=None, checkPassive=True):
        
        GS = 50/MatQ.ZoS    # Normalized input load to the network 
        GL = 50/MatQ.ZoL     # Normalized output load to the network 
        
        # Get MatQ attributes
        M = MatQ.M
        extraNodesS = MatQ.extraNodesS
        extraNodesL = MatQ.extraNodesL
        nodeScalingFactor = MatQ.nodeScalingFactor
            
        # Size of the matrix
        N = np.size(M,0)
        
        # Evaluation points.
        # Since losses are now included in the diag of M, we have to evaluate in the imaginary axis
        if evalS is not None:  
            eval = evalS
        else:  
            eval = 1j * self.s.imag

        # Initialize S-parameter vectors with zeros
        self.S11M = np.zeros(len(eval), dtype='complex128')
        self.S21M = np.zeros(len(eval), dtype='complex128')
        self.S22M = np.zeros(len(eval), dtype='complex128')
        
        condY = 0   # Maximum condition number of matrices to invert
        
        # For each evaluation point find S11 and S21
        for i in range(len(eval)):
            # Computation of Y
            I = np.diag(np.concatenate((np.zeros(extraNodesS), nodeScalingFactor[extraNodesS:-extraNodesL]**2 , np.zeros(extraNodesL))))
            Y = 1j*M + eval[i] * I + np.diag(np.concatenate(([GS], np.zeros(N-2) , [GL])))
            condY = max(condY,  np.linalg.cond(Y)) 
            
            # Computation of Z=inv(Y)
            Z = np.linalg.inv(Y)
            
            # Computation of S11 and S22
            # Sign changed to agree with sign of S11 from CP
            self.S11M[i]= - (1.-2.*GS*Z[0, 0])
            self.S22M[i]= - (1.-2.*GL*Z[-1, -1])
            
            # Computation of S21
            # Sign changed to agree with sign of S11 from CP
            self.S21M[i]=-2.*sqrt(GS*GL)*Z[-1, 0]

        if condY > 10**12:
            warningMsg("WARNING: Computation of [S] parameters from CM:<p>Condition number of matrix to invert is %.1e.<br>[S] parameter values will be wrong due to lack of precision." % condY)
            
        if checkPassive:
            epstol = 1e-3
            S11max,  S21max,  S22max = max(abs(self.S11M)), max(abs(self.S21M)), max(abs(self.S22M))
            if max([S11max,  S21max,  S22max]) > 1+epstol: 
                warningMsg("""WARNING: [S] parameters computed from CM '%s' are larger than 0 dB at some frequency points<p>
                        max[S11] = %.2f dB, max[S21] = %.2f dB, max[S22] = %.2f dB
                      """ % (MatQ.name, 20*log10(S11max), 20*log10(S21max), 20*log10(S22max)) )
        
        ##
        # Phase of S11 return losses computed from coupling matrix
        self.phaseS11M = np.unwrap(np.angle(self.S11M))
        
        ##
        # Phase of S22 return losses computed from coupling matrix
        self.phaseS22M = np.unwrap(np.angle(self.S22M))
        
        ##
        # Phase of transfer function computed from coupling matrix
        self.phaseS21M = np.unwrap(np.angle(self.S21M))
        
        if evalS is None: # Compute group delay only for self.freq frequency axis, which has been defined when computing [S] from polynomials
            ##
            # Group delay computed from coupling matrix 
            self.groupDelayM = -np.diff(self.phaseS21M) / (2*pi* np.diff(self.freq))
    
            # Set causality principle. There may be spurious negative peaks due to computing phase at transmission zeros.
            self.groupDelayM[ self.groupDelayM<0 ] = 0 # CE prefers that the GD goes negative but TAS-E prefers truncation of negative values
            
            ##
            # Frequency axis with sampling points in the middle of self.freq sampling.
            # It is necessary to plot self.groupDelayM
            self.freq2 = (self.freq[:-1] + self.freq[1:])/2
        
            ##
            # Group delay computed from Characteristic Polynomials, sampled at self.freq2
            # It is necessary to compare with self.groupDelayM
            self.groupDelay2 = (self.groupDelay[:-1] + self.groupDelay[1:])/2
            
    
        if evalS is None:  # Computation of different error measures in S11 and S21 parameters
        
            errorS11dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S11)) - 20.*np.log10(np.abs(self.S11M))))
            errorS22dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S22)) - 20.*np.log10(np.abs(self.S22M))))
            errorS21dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S21)) - 20.*np.log10(np.abs(self.S21M))))
            errorGroupDelay = np.mean(np.abs( self.groupDelayM / self.groupDelay2 - 1 ))
        
#        errorS11dBmax = np.max(np.abs(20.*np.log10(np.abs(self.S11))-20.*np.log10(np.abs(self.S11M))))
#        errorS11linavg = np.mean(np.abs(self.S11 - self.S11M)/np.abs(self.S11))
#        errorS11norm = np.linalg.norm(self.S11 - self.S11M) / np.linalg.norm(self.S11)
#        errorS11linmax = np.max(np.abs(self.S11 - self.S11M)/np.abs(self.S11))
#        myPrint("Error S11 (Avg[S11cp(dB)-S11cm(dB)]): %g dB" % errorS11dBavg )
#        myPrint("Error S11 (Avg[|S11cp-S11cm|/|S11cp|]): %g" % errorS11linavg )
#        myPrint("Error S11 ||S11cp-S11cm||/||S11cp||): %g" % errorS11norm )
#        myPrint("Error S11 (Max[S11cp(dB)-S11cm(dB)]): %g dB" % errorS11dBmax )
#        myPrint("Error S11 (Max[|S11cp-S11cm|/|S11cp|]): %g" % errorS11linmax )
            
#        errorS21dBmax = np.max(np.abs(20.*np.log10(np.abs(self.S21))-20.*np.log10(np.abs(self.S21M))))
#        errorS21norm = np.linalg.norm(self.S21 - self.S21M) / np.linalg.norm(self.S21)
#        errorS21linavg = np.mean(np.abs(self.S21 - self.S21M)/np.abs(self.S21))
#        errorS21linmax = np.max(np.abs(self.S21 - self.S21M)/np.abs(self.S21))
#        myPrint("Error S21 (Avg[S21cp(dB)-S21cm(dB)]): %g dB" % errorS21dBavg )
#        myPrint("Error S21 (Avg[|S21cp-S11cm|/|S21cp|]): %g" % errorS21linavg )
#        myPrint("Error S21 ||S21cp-S21cm||/||S21cp||): %g" % errorS21norm )
#        myPrint("Error S21 (Max[S21cp(dB)-S21cm(dB)]): %g dB" % errorS21dBmax )
#        myPrint("Error S21 (Max[|S21cp-S11cm|/|S21cp|]): %g" % errorS21linmax )


#            stErrorCM = """   Average error in [S] computed from Coupling Matrix vs. [S] computed from Characteristic polynomials: 
#            <blockquote>S<sub>11</sub> error = %.3g dB</blockquote>
#            <blockquote>S<sub>22</sub> error = %.3g dB</blockquote>
#            <blockquote>S<sub>21</sub> error = %.3g dB</blockquote>
#            <blockquote>Group delay mean relative error = %.9e</blockquote>""" % (errorS11dBavg,  errorS22dBavg, errorS21dBavg, errorGroupDelay  )
            
            stErrorCM = """   Average error in [S] computed from Coupling Matrix vs. [S] computed from Characteristic polynomials: 
            <blockquote>S<sub>11</sub> error = %.3g dB, S<sub>22</sub> error = %.3g dB, S<sub>21</sub> error = %.3g dB</blockquote>
            <blockquote>Group delay mean relative error = %.3e</blockquote>""" % (errorS11dBavg,  errorS22dBavg, errorS21dBavg, errorGroupDelay  )
        
            maxErrorDB = max([errorS11dBavg,  errorS21dBavg])   # Do not use errorS22dBavg, since errors in S22 may be larger 
            CM_error = abs( 10** (-maxErrorDB/ 20) - 1 )
        
            return stErrorCM, CM_error
        else:
            return None
        
        
    ##
    # 
    #         Compute stored energy and dissipated power at resonators and lossy couplings, from a coupling matrix.
    #         Evaluation points are taken from attribute Sparameters.s .
    #         Output in atributes Sparameters.energyRes, Sparameters.powerRes , 
    #         which are 2D numpy arrays, each column is the energy(f) or power(f) at a resonator.
    #         @param MatQ = MatrixQ class instance, containing the coupling matrix to process.
    #         
    #         Dissipated power at \f$i\f$-th resonator can be found as:
    #         \f[
    #             P_i(f) = \frac{G_i}{2} |V_i(f)|^2 
    #         \f]
    #         being \f$ V_i \f$ the voltage at each resonator and \f$ G_i \f$ the conductance of unloaded resonators
    #         \f[
    #         G_i = \frac{Z_0}{Q_i \mathrm{FBW}}
    #         \f]
    #         where \f$ Q_i \f$ is the quality factor, \f$ \mathrm{FBW} = \Delta f / f_0 \f$ is the filter fractional bandwidth and \f$ Z_0 \f$ is the characteristic impedance of the resonators.
    #         
    #         The stored energy is: 
    #         \f[
    #             W_i(f) = \frac {Q_i P_i(f)}{2 \pi f_0} = \frac{Z_0}{4  \pi f_0 \, \mathrm{FBW}} |V_i(f)|^2 
    #         \f]
    #         
    #         For resistive couplings between resonators \f$ i \f$ and \f$ j \f$, with conductance \f$ G_{ij} = \Im M_{ij} \f$, the voltage drop at the coupling is \f$ V_{ij}=V_i-V_j \f$ and, hence,
    #         the dissipated power becomes:
    #         \f[
    #             P_{ij}(f) = \frac{G_{ij}}{2} \; |V_{ij}(f)|^2
    #         \f]
    #         
    #         Following [Pozar eq. (4.14) and (4.16)], the time harmonic complex power is
    #         \f[ P = \frac{1}{2} V I^* = \frac{1}{2} |V|^2 Y^* = Pd + 2 j \omega (W_m - W_e)
    #         \f]
    #         and, accounting that \f$ Y_{ij} = j M_{ij} \f$, the stored energy at magnetic or electric couplings between resonators \f$ i \f$ and \f$ j \f$, 
    #         with susceptance \f$ B_{ij} = \Im Y_{ij} = \Re M_{ij} \f$, is
    #         \f[
    #             W_{ij}(f) = W_m - W_e = \frac{B_{ij}}{4\omega} \; |V_{ij}(f)|^2
    #         \f]
    #         
    #         Voltages at each resonator for unitary current excitation at source port are obtained as:
    #         \f[ V(s)=Z(s) \, J(s) \f]
    #         where \f$ J(s) \f$ is the input current vector   \f$ J(s) = [I_s \;\: 0 \ldots 0]^T \f$  
    #         
    #         If the input is matched, the input current is \f$ I_s = \sqrt{8 P_0 / Z_0} \f$. We will assume input power \f$ P_0 \f$ equal to 1 Watt.
    #     
    #         The impedance matrix is computed as:
    #         \f[ Z(s)= \left( jM+sI+G  \right)^{-1} \f]
    #         with \f$I \f$ a diagonal matrix with the node scaling factor squared in the diagonal except in the non-resonant nodes positions, where it is zero. 
    #         It means putting zeros in the extraNodesS first positions of the diagonal and in the last extraNodesL positions.
    #         \f[ G = \begin{bmatrix} G_{S} & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & & G_{L}\end{bmatrix} \f]
    #         
    def energyPowerFromCM(self,  MatQ):
        # Threshold to identify non-zero couplings
        eps = MatQ.parent.matrixElementsEps('Energy and power')
        
        GS = 50/MatQ.ZoS    # Normalized input load to the network 
        GL = 50/MatQ.ZoL     # Normalized output load to the network
        
        # Input
        Po = 1.      # Reference power 1W
        Zo = MatQ.Zo    # Characteristic impedance of resonators used to compute the band-pass version of this Coupling Matrix.
        Is = sqrt(Po*8/Zo)  # Source current
        
        # Resonators Q and G
        Qi = MatQ.Q
        G = Zo / (Qi * self.FBW)
        
        # Get MatQ attributes
        M = MatQ.M
        extraNodesS = MatQ.extraNodesS
        extraNodesL = MatQ.extraNodesL
        nodeScalingFactor = MatQ.nodeScalingFactor
            
        # Size of the matrix
        N = np.size(M,0)
        
        # Evaluation points.
        # Since losses are now included in the diag of M, we have to evaluate in the imaginary axis
        eval = 1j * np.imag(self.s)
        Nf = len(eval) # Number of frequency points
        
        # Position and number of resonant nodes
        first = extraNodesS
        last = N - extraNodesL
        Nr = last - first

        # Initialize output vectors with zeros
        energyRes, energyResRev = np.zeros((Nf,  Nr)), np.zeros((Nf,  Nr))
        powerRes, powerResRev = np.zeros((Nf,  Nr)), np.zeros((Nf,  Nr))
        
        # Prepare data for  power dissipated at resistive couplings
        powerCoup = []
        for m in range(N):
            for n in range(m+1, N):
                if n==m: continue
                Gmn = M[m, n].imag
                if np.abs(Gmn) > eps:
                    powerCoup.append( ( (m, n), Gmn, np.zeros(Nf),  np.zeros(Nf) ) )
                    
        # Prepare data for energy stored at magnetic or electric couplings
        energyCoup = []        
        for m in range(N):
            for n in range(m+1, N):
                if n==m: continue
                Bmn = M[m, n].real
                if np.abs(Bmn) > eps:
                    energyCoup.append( ( (m, n), Bmn, np.zeros(Nf), np.zeros(Nf)  ) )
        
        condY = 0   # Maximum condition number of matrices to invert
        
        # For each frequency evaluation point
        for i in range(Nf):
            # Computation of Z
            I = np.diag(np.concatenate((np.zeros(extraNodesS), nodeScalingFactor[extraNodesS:-extraNodesL]**2 , np.zeros(extraNodesL))))
            Z = np.linalg.inv( 1j*M + eval[i] * I + np.diag(np.concatenate(([GS], np.zeros(N-2) , [GL]))) )
            condY = max(condY,  np.linalg.cond(Z)) 
            
            # Vi is the first column of Z times the current, with direct excitation, and the last column with reverse excitation
            V = Z[:, 0]*Is
            Vrev = Z[:, -1]*Is
            
            V2 = np.abs(V[first:last])**2
            Vrev2 = np.abs(Vrev[first:last])**2
            
            powerRes[i, :] = V2 * G/2
            powerResRev[i, :] = Vrev2 * G/2
            
            energyRes[i, :] = V2 * Zo / (4*pi*self.f0*self.FBW)
            energyResRev[i, :] = Vrev2 * Zo / (4*pi*self.f0*self.FBW)
            
            # Compute power dissipated at resistive couplings
            for coup in powerCoup:
                (m, n) = coup[0]
                Gmn = coup[1]
                coup[2][i] = np.abs(V[m] - V[n])**2 * Gmn / 2
                coup[3][i] = np.abs(Vrev[m] - Vrev[n])**2 * Gmn / 2
                
            # Compute energy stored at magnetic or electric couplings
            for coup in energyCoup:
                (m, n) = coup[0]
                Bmn = coup[1]
                coup[2][i] = np.abs(V[m] - V[n])**2 * Bmn / (4* 2*pi*self.freq[i])
                coup[3][i] = np.abs(Vrev[m] - Vrev[n])**2 * Bmn / (4* 2*pi*self.freq[i])
            
        if condY > 10**12:
            warningMsg("""WARNING: Computation of [S] parameters from CM:<p>Condition number of matrix to invert is %.1e.<br>Results will be wrong due to lack of precision.""" % condY)
            
        self.powerRes = powerRes
        self.powerResRev = powerResRev
        self.energyRes = energyRes
        self.energyResRev = energyResRev
        self.powerCoup = powerCoup
        self.energyCoup = energyCoup
        
     
    ##
    # 
    #         Save [S] parameters in output file.
    #         S11, S12, S21 and S22 are attributes of sParameters class, and of type numpy.array.
    #             @param P = Parameter class instance.
    #             @param prec = Number of correct significant digits to print.
    #         
    def saveSparam(self,  P,  prec):
        
        fileName = P.outDirName
        fout = open(fileName + '.s2p', 'w')
        fcsv = open(fileName + '.csv', 'w')
        fd = csv.writer(fcsv, delimiter = ';')
        
        # Set column width
        colWidth = 10 + prec
        
        fout.write('! %s; version %s\n' % pkgNameVersion(__doc__) )
        fout.write ('! [S] parameters file in TouchStone format.\n\n')
        fout.write('# GHZ S DB R 50\n\n')
        
        fd.writerow(pkgNameVersion(__doc__))
        fd.writerow(['[S] parameters'])
        
        st = '%-*s%-*s%-*s%-*s%-*s%-*s%-*s%-*s%-*s\n' % (colWidth, '  !Freq', colWidth, '  S11(dB)', colWidth, '  S11(ang)',  colWidth, '  S12(dB)', colWidth, '  S12(ang)',  colWidth, '  S21(dB)', colWidth, '  S21(ang)',  colWidth, '  S22(dB)',  colWidth, '  S22(ang)' )
        fout.write(st)
        fd.writerow(['Freq(GHz)', 'S11(dB)', 'S11(ang)', 'S12(dB)', 'S12(ang)', 'S21(dB)', 'S21(ang)',  'S22(dB)',  'S22(ang)'])
        
        st = ''
        for n in range(len(self.S11)):
            freq = str(self.freq[n] / 1e9)
            S11 = complexStr(self.S11[n],  prec,  'DB', colWidth)
            S12 = complexStr(self.S12[n],  prec, 'DB',  colWidth)
            S21 = complexStr(self.S21[n],  prec,  'DB', colWidth)
            S22 = complexStr(self.S22[n],  prec,  'DB', colWidth)
            st += '%-*s%-*s%-*s%-*s%-*s\n' % (colWidth, freq, colWidth,  S11,  colWidth,  S21,  colWidth,  S12,  colWidth,  S22)
            fd.writerow([freq] + complexStr(self.S11[n],  prec,  'DB', colWidth, listFormat = True) + complexStr(self.S12[n],  prec,  'DB', colWidth, listFormat = True) + complexStr(self.S21[n],  prec,  'DB', colWidth, listFormat = True) + complexStr(self.S22[n],  prec,  'DB', colWidth, listFormat = True))
        fout.write(st)
        fout.close()
        fcsv.close()
        
        myPrint("[S] parameters saved to files: %s.s2p and %s.csv" % (fileName,  fileName) )
        
        
# End class: Sparameters
#!!!!!!!!!!!!!!!!!!!!!!!!!#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin class: M

##
# 
#     Class that contains a coupling matrix, its name, the number of non-resonant nodes, the Q of resonators and some operations.
#     
class MatrixQ(object):
    
    ##
    # 
    #         Initialize MatrixQ class. Computes Q of resonators and stores the vector of Q values in self.Q .
    #         @param parent = CouplingMatrices class instance to which this coupling matrix belongs.
    #         @param M = Coupling matrix.
    #         @param name = Matrix name.
    #         @param extraNodesS = Number of extra rows and columns in current matrix M at the source side (non-resonant nodes).
    #         @param extraNodesL = Number of extra rows and columns in current matrix M at the load side (non-resonant nodes).
    #         @param FT = Frequency transformation class instance.
    #         @param fileExt = Extension of file name, if this matrix is to be automatically saved in a file (with the output file name in the parameters class). 
    #                           If equal to None, this matrix is not automatically saved. Anyway the matrix can be manually saved. Default None.
    #         @param nodeScalingFactor = Vector containing in each node position the factor of the node scaling applied to that node. Default None.
    #         @param flagRotateAllowed = Flag that indicates if it is allowed to do a matrix rotation (Impossible after a resonant node scaling). Default True.
    #         
    def __init__(self, parent, M,  name,  extraNodesS,  extraNodesL,  FT,  fileExt = None,  nodeScalingFactor=None,  flagRotateAllowed = True):
        ##
        # Coupling matrix
        self.M = M.copy()   # Deep copy of numpy array
        
        ##
        # Matrix order
        self.order = np.size(M, 0) 
        
        ##
        # Number of extra rows and columns in current matrix M at the source side (non-resonant nodes).
        self.extraNodesS = extraNodesS

        ##
        # Number of extra rows and columns in current matrix M at the load side (non-resonant nodes).
        self.extraNodesL = extraNodesL
        
        ##
        # Matrix name.
        # It is converted to unicode to allow safe comparison with matrix names entered by the user in Matrix Edit Gui.
        self.name =  unicode(name)
        
        ##
        # Frequency transform to use when computing Q of resonators.
        self.FT = FT
        
        ##
        # Extension of file name, if this matrix is to be saved in a file. Otherwise None.
        self.fileExt = fileExt
        
        ##
        # CouplingMatrices class instance to which this coupling matrix belongs.
        self.parent = parent
        
        ##
        # Band pass version of the matrix
        self.Mbp = None
        
        ##
        # Coupling bandwidth version of the matrix
        self.Mcbw = None
        
        ##
        # Coupling bandwidth with transmission lines version of the matrix
        self.McbwTL = None
        
        ##
        # Labels for rows and columns. Will be used for saving energy and power.
        self.rowIndices = None
        
        ##
        # Characteristic impedance of resonators
        self.Zo = 50.
        
        ##
        # Characteristic impedance of source
        self.ZoS = 50.
        
        ##
        # Characteristic impedance of load
        self.ZoL = 50.

        # Compute Q of resonators
        self.Qresonators()
       
        # Vector with the information of the performed node scalings.
        if nodeScalingFactor is not None:
            self.nodeScalingFactor = nodeScalingFactor.copy()
        else:
            self.nodeScalingFactor = np.ones(np.size(self.M, 0))
        
        # Flag to control if it is possible to do any matrix rotation.
        self.flagRotateAllowed = flagRotateAllowed

        
    ##
    # 
    #         Function to copy instances of the MatrixQ class (copies the self instance).
    #             @return New instance of MatrixQ class copied from the self instance.
    #         
    def copy(self):
        return MatrixQ(self.parent, self.M,  self.name,  self.extraNodesS,  self.extraNodesL,  self.FT,  fileExt = self.fileExt, nodeScalingFactor= self.nodeScalingFactor, flagRotateAllowed = self.flagRotateAllowed)
            
    
    ##
    # 
    #         Function to set name field of the instance of the MatrixQ class.
    #             @param name = Matrix name.
    #         
    def setName(self, name):
        self.name = name
        
    
    ##
    # 
    #         Function to set file extension field of the instance of the MatrixQ class.
    #             @param fileExt = Extension of file name, if this matrix is to be saved in a file. Otherwise None.
    #         
    def setFileExt(self, fileExt):
        self.fileExt = fileExt
        
    
    ##
    # 
    #         Compute the quality factors of each resonant node of coupling matrix self.M.
    #         Frequency transform is self.FT and output Q is stored in self.Q.
    #         
    #         For each node \f$ n \f$, the quality factor is 
    #         \f[
    #              Q_{n} = \frac{1}{G_n \, \mathrm{FBW}} 
    #          \f]
    #         where the loss conductance of the unloaded resonator is:
    #         \f[
    #             G_{n} = - \sum_{k} \Im M_{kn}
    #         \f]
    #      
    #         <b> Proof: </b>
    #         
    #         In the parallel N+2 normalized lowpass prototype equivalent circuit, at each node we have current excitation \f$ I_n \f$ equal to 
    #         the current that flows away from the node through the circuit elements and couplings:
    #         \f[
    #             I_n = \left( G_n + j B_n + s C_n \right) V_n + \sum_{k \neq n} j M'_{nk} V_k + \sum_{k \neq n} G_{nk} (V_n - V_k)
    #          \f]
    #         where \f$ M'_{kn} \f$ are the real couplings, \f$ G_n \f$ is the unloaded resonator conductance (losses), \f$ G_{kn} \f$ is the conductance of resistive
    #         couplings between resonators and, obviously, the excitation is  \f$ I_n = 0 \f$ for  \f$ n \neq S \f$ . For \f$ n = S \f$ we have \f$ I_n = I_S \f$ (source) .
    #         
    #         If we define the mutual admittances between nodes as:
    #         \f[
    #             Y_{nm} = \left. \frac{I_n}{V_m} \right|_{V_{k \neq m}= 0} 
    #          \f]
    #         we have 
    #         \f[  
    #             Y_{nn} = G_n + j B_n + s C_n + \sum_{k \neq n} G_{kn}
    #          \f]
    #          and
    #          \f[  
    #             Y_{nm} = j M'_{nm} - G_{nm}
    #          \f]
    #         In matrix form, with \f$ C_n = 1 \f$ :
    #         \f[
    #             [Y] = j[M] + s [I]
    #          \f]
    #         where the elements of the coupling matrix \f$ [M] \f$ are:
    #         \f[
    #             M_{nn} = \frac{Y_{nn} - s}{j} = B_n - j G_{n} - j \sum_{k \neq n} G_{kn}
    #          \f]
    #         and 
    #         \f[
    #             M_{nm} = \frac{Y_{nm}}{j} = M'_{nm} + j G_{nm}
    #         \f]
    #         Therefore:
    #         \f[
    #             \Im M_{nn} = - G_{n} - \sum_{k \neq n} G_{kn} = - G_{n} - \sum_{k \neq n} \Im M_{kn}
    #         \f]
    #         and we can compute the resonator unloaded conductance as:
    #         \f[
    #             G_{n} = - \Im M_{nn} - \sum_{k \neq n} \Im M_{kn} = - \sum_{k} \Im M_{kn}
    #         \f]
    #         
    def Qresonators(self):
        eps = self.parent.matrixElementsEps()     # Threshold to determine that losses are zero in presence of roundoff errors
        FBW = self.FT.FBW
        M = self.M
        first = self.extraNodesS
        last = np.size(M, 1) - self.extraNodesL
        
        losses = np.zeros(last - first)
        for i in range(first,  last):
            losses[i-first] = sum(np.imag(M[:,i]))
            if abs(losses[i-first]) < eps:
                losses[i-first] = 0.
        
        ##
        # Q of resonators
        self.Q = -1 / (losses * FBW)
        
        
    ##
    # 
    #         Compute bandpass and coupling bandwidth versions of the matrix.
    #         
    #         <b>Coupling bandwidth matrix</b>
    #         
    #         The diagonal elements of the coupling bandwidth matrix \f$ [\Delta f] \f$ is obtained as:
    #           \f[
    #             \Delta f_{ii} = f_i
    #           \f]
    #           where the unnormalised resonant frequencies (Hz) of the resonators are:
    #           \f[
    #             f_i = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
    #           \f]
    #           with \f$ \omega'' = - B_i \Delta f / f_0 \f$ and \f$ B_i = \Re M_{ii} \f$ . The factor \f$ \Delta f / f_0 \f$ is the fractional bandwidth.
    # 
    #           For non-resonant nodes \f$ f_i = f_0 \f$ .
    #           
    #         The couplings with source and load are:
    #           \f[
    #             \Delta f_{Si} = M_{Si}^2 \Delta f \qquad \Delta f_{Li} = M_{Li}^2 \Delta f
    #           \f]
    #           \f[
    #             \Delta f_{SS} =  \Delta f_{LL} = 0
    #           \f]
    #           where \f$ M_{ij} \f$ are the normalised low-pass coupling matrix elements and \f$ \Delta f \f$ is the bandwidth (MHz).
    #           
    #           The couplings between inner nodes (resonators or non-resonant nodes) are:
    #           \f[
    #             \Delta f_{ij} = M_{ij} \Delta f
    #           \f]
    #         
    #         <b>Coupling bandwidth matrix with transmission lines </b>
    #      
    #         For each matrix element \f$ i, j \f$, the scaling factor for the pass-band coupling matrix is 
    #         \f[
    #             K_{ij} = \alpha_{ij} M_{ij}
    #         \f]
    #         where \f$ [M] \f$ is the low-pass normalised matrix and \f$ [K] \f$ is the pass-band coupling matrix.
    #         \f[
    #              \alpha_{ij} = \frac{\pi}{2} \frac{\mathrm{FBW}}{Z_0} \qquad \mathrm{if} \; i \; \mathrm{and} \; j \neq S \; \mathrm{or} \; L
    #          \f]
    #          \f[
    #              \alpha_{ij} = \sqrt{\frac{\pi}{2} \frac{\mathrm{FBW}}{Z_0 Z_{S,L}}} \qquad \mathrm{if} \; i \; \mathrm{or} \; j = S \; \mathrm{or} \; L
    #          \f]   
    #          \f[
    #             \alpha_{SS} = \alpha_{SL} = \alpha_{LS} = \alpha_{LL} = 1
    #          \f]
    #          \f[
    #              Z_0 = Z_S = Z_L = 50  \Omega
    #           \f]
    #           
    #           The diagonal elements of the coupling bandwidth matrix \f$ [\Delta f] \f$ is obtained as:
    #           \f[
    #             \Delta f_{ii} = f_i
    #           \f]
    #           where the unnormalised resonant frequencies (Hz) of the resonators are:
    #           \f[
    #             f_i = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
    #           \f]
    #           with \f$ \omega'' = - B_i \Delta f / f_0 \f$ and \f$ B_i = \Re M_{ii} \f$ . The factor \f$ \Delta f / f_0 \f$ is the fractional bandwidth.
    # 
    #           For non-resonant nodes \f$ f_i = f_0 \f$ .
    #           
    #           The couplings with source and load are:
    #           \f[
    #             \Delta f_{Si} = \frac{f0}{K_{Si}} \qquad \Delta f_{Li} = \frac{f0}{K_{Li}}
    #           \f]
    #           \f[
    #             \Delta f_{SS} =  \Delta f_{LL} = 0
    #           \f]
    #           
    #           The couplings between inner nodes (resonators or non-resonant nodes) are:
    #           \f[
    #             \Delta f_{ij} = f_0 \sqrt{ \left( \frac{K_{ij}}{A} \right) ^2 + C }
    #           \f]
    #           with
    #           \f[
    #             A = \frac{1}{2} \left( \frac{f_j}{f_i} + \frac{f_i}{f_j} \right)
    #           \f]
    #           \f[
    #             C = \left| \frac{f_j^2 - f_i^2}{f_j^2 + f_i^2} \right| ^2
    #           \f]
    # 
    #           For syncronous nodes, \f$ f_i = f_0 \f$ and therefore
    #           \f[
    #             A = 1 \qquad C = 0
    #           \f]therefore
    #           \f[
    #             \Delta f_{ij} = K_{ij} f_0
    #           \f]
    #         
    def lp2bp2cbw(self):
        eps = self.parent.matrixElementsEps()     # Threshold to determine that couplings are zero in presence of roundoff errors
        
        N = np.size(self.M, 1)
        FBW = self.FT.FBW
        BW = self.FT.BW
        f0 = self.FT.f0
        
        alpha = np.matrix(np.ones(N))
        
        # Jordi Mateu's approach, with transmission lines
        alpha[0,  1:N-1] = sqrt( self.FT.FBW*pi/(2*self.Zo) )   # alpha is a 2-D matrix
        alpha[0,  0] = sqrt(1./self.ZoS)     # ZoS is the source impedance, normally 50 ohm
        alpha[0,  N-1] = sqrt(1./self.ZoL)     # ZoL is the source impedance, normally 50 ohm
        self.Mbp = self.M * np.array((alpha.T*alpha))
        
        # Just M*FBW, as in CMSDemo software
#        self.Mbp = self.M * FBW
        
        self.Mcbw = np.zeros((N, N), dtype=complex)
        self.McbwTL = np.zeros((N, N), dtype=complex)
        
        # Resonant frequencies in Hz:
        Bn_bp = np.diag(self.M).real * FBW
        fn = self.FT.unormFreq(-np.diag(self.M).real)
        fn2 = fn*fn
        
        for m in range(N):
            for n in range(m+1):
                if abs(self.M[m, n]) < eps: continue    # Null coupling
                if m==N-1 or n==0:  # "L" row or "S" column, and their symmetrical.
                    self.McbwTL[m, n] = (f0/1e6) / self.Mbp[m, n]
                    self.McbwTL[n, m] = self.McbwTL[m, n]
                    self.Mcbw[m, n] = self.M[m, n]**2 * BW / 1e6
                    self.Mcbw[n, m] = self.Mcbw[m, n]
                else:   # Couplings between "R" and/or "NR" nodes
                    if n==m:    # Main diagonal
                        self.McbwTL[n, n] = fn[n] / 1e6   # MHz
                        self.Mcbw[n, n]    = fn[n] / 1e6   # MHz
                    else:         # Couplings
                        A = ( fn[n]/fn[m] + fn[m]/fn[n] ) / 2
                        C = abs( ( (fn2[n]-fn2[m])/(fn2[m]+fn2[m]) )**2 )
#                        print '(m,n) =', (m, n)
#                        print 'Kij =', self.Mbp[m, n]
#                        print 'A, C, (Kij/A)**2, sqrt((Kij/A)**2 + C) =', A,  C,  (self.Mbp[m, n]/A)**2, np.sqrt( (self.Mbp[m, n]/A)**2 + abs(C) )
                        self.McbwTL[m, n] = f0 * np.sqrt( (self.Mbp[m, n]/A)**2 + C ) / 1e6   # MHz
                    
                        # Check if we have to take the negative sign of the sqrt, to preserve the same sign as in self.Mbp[m, n]*f0
                        if abs(self.McbwTL[m, n] / (self.Mbp[m, n]*f0/1e6) + 1) < 1: self.McbwTL[m, n] = - self.McbwTL[m, n]
                        
                        self.McbwTL[n, m] = self.McbwTL[m, n]
                        
                        self.Mcbw[m, n] = self.M[m, n] * BW / 1e6
                        self.Mcbw[n, m] = self.Mcbw[m, n]
    
    
    ##
    # 
    #         Rotate the coupling matrix self.M. The resulting self.M matrix maintains the eigenvalues of the original matrix and his reflection and transmission properties.
    #         @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
    #         @param i = Row index of the pivot that will be used in the rotation.
    #         @param j = Column index of the pivot that will be used in the rotation.
    #         @param rotAng = Rotation angle (rad).
    #         @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
    #         
    #         [TN 102.2 sec. 9.1]
    #         By definition, a rotation matrix is a matrix where every value in the diagonal is 1 except the two values \f$R_{ii}=R_{jj}=c_{r}\f$ and the rest of the values of the matrix are zero except the two values \f$R_{ji}=-R_{ij}=s_{r}\f$, where  [TN 102.2 eq.(A.2) and (A.3)]: 
    #         \f[\begin{matrix}[i,j] & = & \text{rotation pivot} &\\ 
    #             c_r & = & cos(rotAng), & \text{if rotationType = 'trigonometric'} \\
    #                  &    &  cosh(rotAng), & \text{if rotationType = 'hyperbolic'} \\
    #             s_r & = & sin(rotAng), & \text{if rotationType = 'trigonometric'} \\
    #                  &    &  jsinh(rotAng), & \text{if rotationType = 'hyperbolic'} \end{matrix}
    #         \f]
    #         
    #         To apply the rotation we just need to do the product [TN 102.2 eq.(A.1)]:
    #         \f[M_{rot} = RMR^{T}\f]
    #         where \f$R\f$ is the rotation matrix described above.
    #         
    #         
    def rotateMatrix(self, rotationType, i, j, rotAng, flagQ=True):
        
        M = np.array(self.M)
        N = np.size(M,0)        # Matrix dimension
        
        if not self.flagRotateAllowed:
            raise synthError,  u"It is not allowed to rotate matrix after a Node Scaling in a resonant node with scale factor different from 1 or -1."
            
        # Row and column header labels
        labels = ['S']
        for n in range(0, self.extraNodesS-1 ): labels.append('NR')
        labels.extend(  [ 'R' for elem in range(0,  N - self.extraNodesS - self.extraNodesL) ] )
        for n in range(0,  self.extraNodesL-1 ): labels.append('NR')
        labels.append('L')
            
        # Check cases of forbidden rotations
        if i == j:
            raise synthError,  u"It is not allowed to rotate a matrix with a pivot in the diagonal."
        elif labels[i] == 'S' or labels[j] == 'S' or labels[i] == 'L' or labels[j] == 'L':
            raise synthError,  u"It is not allowed to rotate a matrix with a pivot in the S or L rows or columns."
        elif labels[i]  != labels[j]:
            raise synthError,  u"It is not allowed to rotate a matrix with a pivot in a row and column of different resonant/non-resonant type." 
        
        # Computation of the rotation angle
        if rotationType == 'trigonometric':
            cr = np.cos(rotAng)
            sr = np.sin(rotAng)
        elif rotationType == 'hyperbolic':
            cr = np.cosh(rotAng)
            sr = 1j*np.sinh(rotAng)
        else:
            raise synthError,  u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
        
        # Construction of the rotation matrix
        R = np.diag(np.ones(N)+1j*np.zeros(N))
        R[i,i] = cr
        R[j,j] = cr
        R[i,j] = -sr
        R[j,i] = sr
        
        # Application of the rotation matrix R to the original matrix M
        rotatedM = np.dot(R, np.dot(M, R.T))
        self.M = rotatedM
        if flagQ:
            self.Qresonators()        
        
    
#    def mpRotateMatrix(self,  rotationType, i, j, rotAng, flagQ=True):
#        """
#        Rotate the coupling matrix self.M, using multiprecision. The resulting self.M matrix maintains the eigenvalues of the original matrix and his reflection and transmission properties.
#        @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
#        @param i = Row index of the pivot that will be used in the rotation.
#        @param j = Column index of the pivot that will be used in the rotation.
#        @param rotAng = Rotation angle (rad).
#        @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
#        
#        [TN 102.2 sec. 9.1]
#        By definition, a rotation matrix is a matrix where every value in the diagonal is 1 except the two values \f$R_{ii}=R_{jj}=c_{r}\f$ and the rest of the values of the matrix are zero except the two values \f$R_{ji}=-R_{ij}=s_{r}\f$, where  [TN 102.2 eq.(A.2) and (A.3)]: 
#        \f[\begin{matrix}[i,j] & = & \text{rotation pivot} &\\ 
#            c_r & = & cos(rotAng), & \text{if rotationType = 'trigonometric'} \\
#                 &    &  cosh(rotAng), & \text{if rotationType = 'hyperbolic'} \\
#            s_r & = & sin(rotAng), & \text{if rotationType = 'trigonometric'} \\
#                 &    &  jsinh(rotAng), & \text{if rotationType = 'hyperbolic'} \end{matrix}
#        \f]
#        
#        To apply the rotation we just need to do the product [TN 102.2 eq.(A.1)]:
#        \f[M_{rot} = RMR^{T}\f]
#        where \f$R\f$ is the rotation matrix described above.
#        
#        """
#
#        if not self.flagRotateAllowed:
#            raise synthError,  u"It is not allowed to rotate matrix after a Node Scaling in a resonant node."
#        
#        M = np2mp(self.M)
#        
#        # Matrix dimension
#        N = np.size(self.M,0)
#        
#        # Computation of the rotation angle
#        if rotationType == 'trigonometric':
#            cr = mp.cos(rotAng)
#            sr = mp.sin(rotAng)
#        elif rotationType == 'hyperbolic':
#            cr = mp.cosh(rotAng)
#            sr = 1j*mp.sinh(rotAng)
#        else:
#            raise synthError,  u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
#        
#        # Construction of the rotation matrix
#        R = mp.eye(N)
#        R[i,i] = cr
#        R[j,j] = cr
#        R[i,j] = -sr
#        R[j,i] = sr
#        
#        # Application of the rotation matrix R to the original matrix M
#        rotatedM = R * M * R.T
#        self.M = mp2np(rotatedM)
#        
#        if flagQ:
#            self.Qresonators()
    
    
    ##
    # 
    #         Add external non-resonant nodes to the coupling matrix self.M at the source or load sides, or both.
    #         The resulting self.M matrix will have one or two columns and rows more than the input one.
    #         @param position = Position where we want to inflate the matrix (string). It can be 'source', 'load' or 'both'.
    #         
    #         If position == 'source'
    #         \f[
    #         Minflate=\begin{bmatrix} 0 & k_1 & 0 & \hdots  \\ k_1 &  &  &   \\ 0 &  & M &   \\ \vdots &  &  &    \end{bmatrix}
    #         \f]
    #         
    #         If position == 'load'
    #         \f[
    #         Minflate=\begin{bmatrix}    &  &  & \vdots \\   & M &  & 0 \\   &  &  & k_2 \\  \hdots & 0 & k_2 & 0 \end{bmatrix}
    #         \f]
    #         
    #         If position == 'both'
    #         \f[
    #         Minflate=\begin{bmatrix} 0 & k_1 & 0 & \hdots & 0 \\ k_1 &  &  &  & \vdots \\ 0 &  & M &  & 0 \\ \vdots &  &  &  & k_2 \\ 0 & \hdots & 0 & k_2 & 0 \end{bmatrix}
    #         \f]
    #         
    #         where \f$k_1\f$ and \f$k_2\f$ are usually equal to 1 unless the first or last nodes of the matrix have been previously scaled by respective values \f$k_1\f$ and \f$k_2\f$ with node-scaling operations. 
    #         If this is the case, the value of \f$k_1\f$ and \f$k_2\f$ for that nodes are retrieved from the self.nodeScalingFactor attribute.
    #         
    def inflateMatrix(self, position):
        M = self.M
        # Order of the matrix
        N = np.size(M, 0)
        
        if position ==  u'source':
            Minflate = np.zeros((N+1, N+1), dtype='complex128')
            Minflate[1::, 1::] = M
            k = self.nodeScalingFactor[0]
            Minflate[0, 1] = Minflate[1, 0] = k
            self.extraNodesS = self.extraNodesS + 1
            self.nodeScalingFactor = np.concatenate(([1], self.nodeScalingFactor))
        elif position ==  u'load':
            Minflate = np.zeros((N+1, N+1), dtype='complex128')
            Minflate[:-1:, :-1:] = M
            k = self.nodeScalingFactor[-1]
            Minflate[-1, -2] = Minflate[-2, -1] = k
            self.extraNodesL = self.extraNodesL + 1
            self.nodeScalingFactor = np.concatenate(( self.nodeScalingFactor, [1]))
        elif position ==  u'both':
            Minflate = np.zeros((N+2, N+2), dtype='complex128')
            Minflate[1:-1, 1:-1] = M
            k = self.nodeScalingFactor[0]
            Minflate[0, 1] = Minflate[1, 0] = k
            k = self.nodeScalingFactor[-1]
            Minflate[-1, -2] = Minflate[-2, -1] = k
            self.extraNodesS = self.extraNodesS + 1
            self.extraNodesL = self.extraNodesL + 1
            self.nodeScalingFactor = np.concatenate(([1],  self.nodeScalingFactor, [1]))
        else: # Don't do anything
            Minflate = M
        
        self.M = Minflate
        
    
    ##
    # 
    #         Scale a node of the coupling matrix self.M, by multiplying row [i] and column [i] by a factor. Matrix rotations are not allowed after node scaling of a resonant node.
    #         @param i = Index of row and column which we want to scale.
    #         @param b = Factor to use to scale the ith row and column.
    #         @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
    #         
    def scaleNode(self,  i,  b, flagQ=True):
        self.M[:, i] = self.M[:, i]*b
        self.M[i, :] = self.M[i, :]*b
        
        self.nodeScalingFactor[i] = self.nodeScalingFactor[i]*b
        if (i >= self.extraNodesS) and (i < np.size(self.M,0) - self.extraNodesL): # If it is a resonant node
            if abs(b) != 1: # If the scaling factor is not equal to 1 or -1
                self.flagRotateAllowed = False
        if flagQ:
            self.Qresonators()
    

    ##
    # 
    #         Compute the rotation angle such that, after applying a rotation of type rotationType, the element \f$ M[m,n] \f$ of the matrix self.M is equal to zero.
    #         The element to set to zero [m,n] must be in the row or column of the rotation pivot [i,j] or its transpose [j,i].
    #         @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
    #         @param i = Row index of the pivot that will be used in the rotation.
    #         @param j = Column index of the pivot that will be used in the rotation.
    #         @param m = Row index of the element that will become zero after rotation.
    #         @param n = Column index of the element that will become zero after rotation.
    #         @return rotAng = Rotation angle to remove the element \f$ M[m,n] \f$.
    #         
    #         First of all, and following the rules in the next table [TN 102.2 tab. A.1], we need to compute a value \f$ val \f$ which depends on the position \f$ [m,n] \f$ of the element to eliminate and 
    #         on the pivot \f$ [i,j] \f$:
    #         \f[
    #             \begin{tabular}{|l|l|l|l|l|l|}
    #             \hline
    #              & Element to eliminate & $i$ & $j$ & $k$ &  $val$   \\
    #             \hline
    #             1 & $M_{mn}$ & $m$ & $\neq m,n$ & $n$ & $\frac{M_{ik}}{M_{jk}}$\\
    #             \hline
    #             2 & $M_{mn}$ & $\neq m,n$ & $m$ & $n$ & $-\frac{M_{jk}}{M_{ik}}$\\
    #             \hline
    #             3 & $M_{mn}$ & $n$ & $\neq m,n$ & $m$ & $\frac{M_{ki}}{M_{kj}}$\\
    #             \hline
    #             4 & $M_{mn}$ & $\neq m,n$ & $n$ & $m$ & $-\frac{M_{kj}}{M_{ki}}$\\
    #             \hline
    #             5 & $M_{mm}$ & $m$ & $\neq m$ & - & $\frac{-M_{ij}+\sqrt{M_{ij}^{2}-M_{ii}M_{jj}}}{M_{jj}}$\\
    #             \hline
    #             6 & $M_{mm}$ & $\neq m$ & $m$ & - & $\frac{M_{ij}+\sqrt{M_{ij}^{2}-M_{ii}M_{jj}}}{M_{ii}}$\\
    #             \hline
    #             7 & $M_{mn}$ & $m$ & $n$ & - & $\frac{2M_{ij}}{M_{jj}-M_{ii}}$\\
    #             \hline
    #             8 & $M_{mn}$ & $n$ & $m$ & - & $\frac{2M_{ij}}{M_{jj}-M_{ii}}$\\
    #             \hline
    #             \end{tabular}
    #         \f]
    #         We define \f$f_{l}=1\f$ for the cases {1,2,3,4,5,6} and \f$f_{l}=0.5\f$ for cases {7,8}.
    #         Then, the rotation angle which sets the element \f$[m,n]\f$ to zero after applying the rotation is:
    #         \f[\begin{matrix}
    #             rotAng & = & f_{l}tan^{-1}(val), & \text{if rotationType = 'trigonometric'} \\
    #             rotAng & = & -jf_{l}tan^{-1}(val), & \text{if rotationType = 'hyperbolic'} \\
    #             \end{matrix}
    #         \f]
    #         
    def rotAngleEliminate(self, rotationType, i, j, m, n):
        f=1
        M = self.M
        
        # Val computation
        if (i<>m) & (i<>n) & (j<>m) & (j<>n):
            # If the value to set to zero [m,n] is not in column i, in row i, in column j or in row j, then it cannot be removed
            raise synthError,  u"The element to set to zero [m,n] must be in the row or column of the rotation pivot [i,j] or its transpose [j,i]."
        elif (m<>n) & (i==m) & (j<>m) & (j<>n):
            # Case 1
            k=n 
            val = M[i,k]/M[j,k]
        elif (m<>n) & (i<>m) & (i<>n) & (j==m):
            # Case 2
            k=n
            val = -M[j,k]/M[i,k]
        elif (m<>n) & (i==n) & (j<>m) & (j<>n):
            # Case 3
            k=m
            val = M[k,i]/M[k,j]
        elif (m<>n) & (i<>m) & (i<>n) & (j==n):
            # Case 4
            k=m
            val = -M[k,j]/M[k,i]
        elif (m==n) & (i==m) & (j<>m):
            # Case 5
            if M[j,j]:
                val = (-M[i,j]+np.sqrt(M[i,j]**2-M[i,i]*M[j,j]))/M[j,j]
            else:
                raise synthError,  u"M[j,j] is equal to zero and therefore this pivot should be changed with another one which is different from zero."
        elif (m==n) & (i<>m) & (j==m):
            # Case 6
            if M[i,i]:
                val = (M[i,j]+np.sqrt(M[i,j]**2-M[i,i]*M[j,j]))/M[i,i]
            else:
                raise synthError,  u"M[i,i] is equal to zero and therefore this pivot should be changed with another one which is different from zero."
        else:
            # Cases 7 and 8
            f = 0.5
            val = 2.*M[i,j]/(M[j,j]-M[i,i])
        
        # Rotation angle computation
        if rotationType == 'trigonometric':
            rotAng = f*np.arctan(val)
        elif rotationType == 'hyperbolic':
            # arctanh(val/j) = -j arctan(val)
            rotAng = -1j*f*np.arctan(val)
        else:
            raise synthError,  u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
    
        return rotAng
        
        
    ##
    # 
    #         Add losses to resonators by substracting a factor 1/(Qeff*FBW) to the diagonal of the coupling matrix.
    #         The typical application is to set a prescribed flatness in the ideal response of a lossy filter synthesis
    #         by emulating the non-flat response of a lossy filter of quality factor equal to Qeff.
    #         Qeff = Quality factor of the resonators to emulate.
    #         
    def addLossesQeff(self,  Qeff):
        Geff = 1/(Qeff*self.FT.FBW)
        M = self.M   # Pointer to matrix, to simply code
        N = np.size(M, 0)
        first = self.extraNodesS
        last = N - self.extraNodesL
        
        for n in range(first,  last):
            M[n, n] = M[n, n].real + 1j * (M[n, n].imag - Geff)
            
            
    ##
    # 
    #         Prepare optimization parameters depending on matrix topology and forced signs.
    #         
    #         @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
    #         @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part. 
    #         @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
    #         @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
    #         
    #         @return x0 = Initial guess. Double numpy vector with the values of the matrix entries which are different from zero. 
    #                                    It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
    #         @return NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
    #         @return NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix. 
    #                               The entries corresponding to resonators are False
    #         @return numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
    #         @return bounds = List of tuples (Min,Max) with bounds for all variables to optimize. The maximum abs of any matrix element is set to 4 times the largest coupling, with a minimum of 2. 
    #                          This bound also affects the real part of diagonal elements (normalized resonant frequencies).
    #         
    def setupOptimTopology(self, topologyReal, topologyImag, topologySignReal, topologySignImag):
        
        M = self.M   # Pointer to matrix, to simplify code
        N = np.size(M, 0)
        
        # The matrix must be symmetric, so only the upper triangular part is considered
#        NonZeroReal, NonZeroImag  = map(bool, np.triu(topologyReal)), map(bool, np.triu(topologyImag))  # The result of np.triu is an array of int64, we need an array of Bool
        UpperBool = np.fromfunction(lambda i,j: i<=j, (N, N))
        NonZeroReal,NonZeroImag  = topologyReal*UpperBool, topologyImag*UpperBool
        
        # For the resonating nodes in the diagonal, the only degree of freedom is the real part (with any sign), since the imaginary part will be fixed by the Q of resonators.
        Rindices = range(self.extraNodesS, N-self.extraNodesL)
        NonZeroReal[Rindices, Rindices] = True
        NonZeroImag[Rindices, Rindices] = False
        topologySignReal[Rindices, Rindices] = 0
        
        # For the SS, LL and NR diagonal elements, the only degree of freedom is the imaginary part and must be negative
        NRindices = range(0, self.extraNodesS) + range(-self.extraNodesL, 0)
        NonZeroReal[NRindices, NRindices] = False
        NonZeroImag[NRindices, NRindices] = True
        topologySignImag[NRindices, NRindices] = -1

        # Prepare array of forced signs for non-zero elements, to be used for optimization bounds
        signs = np.concatenate([ topologySignReal[NonZeroReal], topologySignImag[NonZeroImag] ])
        
        # Number of non-zero real couplings, to access the 1st non-zero imag coupling after concatenation of the non-zero real and imaginary vectors
        numReal = len(M[NonZeroReal])

        # Maximum abs of any matrix element is 4 times the largest coupling.
        # Be careful! With this implementation, this bound also affects the real part of diagonal elements (normalized resonant frequencies).
        maxAbsElem = 4 * np.max(np.abs(M))
        if maxAbsElem < 2: maxAbsElem = 2   # A reasonable minimum value is 2
        
        bounds = [] # Optimization boundaries
        for n in range(len(signs)):
            lowerBound, upperBound = -maxAbsElem, maxAbsElem    # Default case
            if signs[n] >0: lowerBound = 0  # Must be positive
            if signs[n] <0: upperBound = 0  # Must be negative
            bounds.append( (lowerBound, upperBound) )  # Imag elements
            
        x0 = np.concatenate([np.real(M[NonZeroReal]), np.imag(M[NonZeroImag])])  # Optimization initial guess
        
#        print 'NonZeroReal:'
#        print NonZeroReal
#        print 'topologySignReal:'
#        print topologySignReal
#        print 'NonZeroImag:' 
#        print NonZeroImag
#        print 'topologySignImag:'
#        print topologySignImag
#        print 'x0, signs, Bounds  =' , zip(x0, signs, bounds)
        
        return x0, NonZeroReal, NonZeroImag, numReal, bounds
        
        
    ##
    # 
    #         Return a new MatrixQ object that is a copy of self, except for the M matrix which is created from the NonZeroValues, NonZeroReal, NonZeroImag and Q parameters.
    #         
    #         @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero. 
    #                                             It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
    #         @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
    #         @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix. 
    #                              The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
    #         @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
    #         @param Q =  Prescribed uniform Q for the resonators.
    #         
    #         @return newMatQ = New MatrixQ object
    #         
    def newMatQFromNonZero(self, NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q):
        
        # Create an empty matrix and fill it with the non-zero values
        N = self.order
        M = np.zeros((N, N), dtype = complex)
        M[NonZeroReal] = NonZeroValues[0:numReal]
        M[NonZeroImag] = M[NonZeroImag] + 1j*NonZeroValues[numReal::]
        M = M + np.triu(M, 1).T     # Symmetrize the matrix
        
        # We impose the Q of the resonators by fixing the imaginary part of the diagonal
        losses = -1/(Q*self.FT.FBW)
        Rindices = range(self.extraNodesS, N-self.extraNodesL)                          # Indices of resonant nodes
        for i in Rindices:  M[i, i] = M[i, i] + 1j * (losses - np.sum(M[:, i].imag))   # Since now M[i,i].imag=0, then np.sum(M[:, i].imag) is equal to the sum of non-diag elements
        
        return MatrixQ(self.parent, M,  self.name,  self.extraNodesS,  self.extraNodesL,  self.FT,  fileExt = self.fileExt, nodeScalingFactor= self.nodeScalingFactor, flagRotateAllowed = self.flagRotateAllowed)
    
    
    ##
    # 
    #         Optimize the coupling matrix to achieve a prescribed uniform Q and [S] parameters as close as possible to those computed from the characteristic polynomials.
    #         The non-zero coupling matrix entries to optimize are specified by topologyReal and topologyImag parameters, 
    #         which are set by the matrix topology widget in the GUI.
    #         
    #         @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #         @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
    #         @param weights = Optimization weights for the difference in S11 and S21 between the computation from the coupling matrix and from polynomials: 
    #                          List [ weightS11,  weightS22, weightS21, weightGD ]
    #         @param Q = Prescribed uniform Q for the resonators.
    #         @param algor = Constrained optimization algorithm. String equal to 'L-BFG-S', 'TNC' or 'SLSQP'. See below.
    #         @param maxIter = Maximum number of iterations, or function evaluations, depending on algorithm.
    #         @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between P.minFreq and P.maxFreq.
    #         @param S11DR = Dynamic range of S11 and S22 to test.
    #         @param S21DR  = Dynamic range of S21 to test.
    #         @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
    #         @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part.
    #         @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
    #         @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
    #         
    #         \section algors Constrained optimization algorithms:
    #         
    #         - <a href="http://en.wikipedia.org/wiki/Optimization"> http://en.wikipedia.org/wiki/Optimization</a>
    #         
    #         - <a href="http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html"> ttp://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html</a>
    #         
    #         In optimization, quasi-Newton methods (also known as variable metric methods) are algorithms for finding local maxima and minima of functions. 
    #         Quasi-Newton methods are based on Newton's method to find the stationary point of a function, where the gradient is 0. 
    #         Newton's method assumes that the function can be locally approximated as a quadratic in the region around the optimum, 
    #         and use the first and second derivatives (gradient and Hessian) to find the stationary point.
    #         In Quasi-Newton methods the Hessian matrix of second derivatives of the function to be minimized does not need to be computed. 
    #         The Hessian is updated by analyzing successive gradient vectors instead (or approximate gradient evaluations computed using finite differences). 
    #         Quasi-Newton methods are a generalization of the secant method to find the root of the first derivative for multidimensional problems. 
    #         In multi-dimensions the secant equation is under-determined (has not a unique solution), 
    #         and quasi-Newton methods differ in how they constrain the solution, typically by adding a simple low-rank update to the current estimate of the Hessian.
    #         
    #         A necessary condition for optimality is that the gradient be zero. 
    #         Quasi-Newton's method need not converge unless the function has a quadratic Taylor expansion near an optimum.
    #         
    #         - <a href="http://en.wikipedia.org/wiki/Quasi-Newton_methods"> http://en.wikipedia.org/wiki/Quasi-Newton_methods</a>
    #         
    #         - <a href="http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization"> http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization</a>
    #         
    #         \subsection lbfgsb Limited memory variation Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS-B)
    #         
    #         Limited memory variation of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update to approximate the inverse Hessian matrix.
    #         
    #         The BFGS method is one of the most popular quasi-Newton optimization methods. L-BFGS is a limited-memory version of BFGS that is particularly suited to problems with very large numbers of variables. 
    #         
    #         References:
    #         
    #         - <a href="http://en.wikipedia.org/wiki/L-BFGS"> http://en.wikipedia.org/wiki/L-BFGS</a>
    #         
    #         - <a href="http://en.wikipedia.org/wiki/BFGS_method"> http://en.wikipedia.org/wiki/BFGS_method</a>
    #         
    #         - Wright, and Nocedal ‘Numerical Optimization’, 1999, pg. 198.
    #         
    #         - R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208.
    #         
    #         - C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550 - 560.
    #         
    #         \subsection slsqp Sequential Least SQuares Programming (SLSQP)
    #         
    #         SLSQP optimizer is a sequential least squares programming algorithm which uses the Han–Powell quasi–Newton method with a BFGS update of the B–matrix and 
    #         an L1–test function in the step–length algorithm. The optimizer uses a slightly modified version of Lawson and Hanson’s NNLS nonlinear least-squares solver.
    #         
    #         Sequential quadratic programming (SQP) is one of the most popular and robust algorithms for nonlinear continuous optimization. 
    #         The method is based on solving a series of subproblems designed to minimize a quadratic model of the objective subject to a linearization of the constraints. 
    #         If the problem is unconstrained, then the method reduces to Newton's method for finding a point where the gradient of the objective vanishes. 
    #         If the problem has only equality constraints, then the method is equivalent to applying Newton's method to the first-order optimality conditions, or Karush–Kuhn–Tucker conditions, of the problem. 
    #         A number of packages (including NPSOL, NLPQL, OPSYC, OPTIMA, MATLAB, and SQP) are founded on this approach.
    #      
    #         References:
    #         
    #         - <a href="http://en.wikipedia.org/wiki/Sequential_quadratic_programming">  http://en.wikipedia.org/wiki/Sequential_quadratic_programming</a>
    #         
    #         - Kraft D (1988) A software package for sequential quadratic programming. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center — Institute for Flight Mechanics, Koln, Germany.   
    #         
    #         \subsection tnc Bound-constrained truncated newton algorithm (TNC)
    #         
    #         Bound-constrained truncated newton algorithm using gradient information. 
    #         
    #         Truncated-Newton methods are a family of methods suitable for solving large nonlinear optimization problems. 
    #         At each iteration, the current estimate of the solution is updated (i.e., a step is computed) by approximately solving the Newton equations using an iterative algorithm. 
    #         This results in a doubly iterative method: an outer iteration for the nonlinear optimization problem, and an inner iteration for the Newton equations. 
    #         The inner iteration is typically stopped or truncated before the solution to the Newton equations is obtained.
    #         Over the past two decades, a solid convergence theory has been derived for the methods. 
    #         In addition, many algorithmic enhancements have been developed and studied, resulting in a number of publicly-available software packages.  
    #         The result has been a collection of powerful, flexible, and adaptable tools for large-scale nonlinear optimization.
    #         
    #         References:
    #         
    #         - <a href="http://js2007.free.fr/code/index.html#TNC"> http://js2007.free.fr/code/index.html#TNC</a>
    #         
    #         - <a href="http://iris.gmu.edu/~snash/nash/assets/TN_Survey/tn_survey.html"> http://iris.gmu.edu/~snash/nash/assets/TN_Survey/tn_survey.html</a>
    #         
    #         - Stephen G. Nash, “A Survey of Truncated-Newton Methods”, Journal of Computational and Applied Mathematics, 124 (2000), pp. 45-59. 
    #         <a href="http://iris.gmu.edu/~snash/papers/nash_3_99.ps">  http://iris.gmu.edu/~snash/papers/nash_3_99.ps</a>
    #         
    #         
    def uniformQ(self, P, CP, weights, Q, algor, maxIter, Nsamples, S11DR, S21DR, topologyReal, topologyImag, topologySignReal, topologySignImag):
        
        #!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        # BEGIN: Helper function
        ##
        # 
        #             Return a SParameters class instance with sampling frequencies to optimize. 
        #             S11, S22, S21 and groupDelay attributes are initialized with those computed from the polynomials.
        #             
        #             @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
        #             @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
        #             @param FT = FrequencyTransform class instance, containing f0 and BW parameters.
        #             @param minFreq = Minimum frequency to optimize.
        #             @param maxFreq = Maximum frequency to optimize.
        #             @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between minFreq and maxFreq.
        #             
        #             @return SPref = SParameters class instance with sampling frequencies to optimize and [S] parameters computed from CP.
        #             
        def setupOptimFreq(P,  CP, FT,  minFreq, maxFreq, Nsamples):
            
            # Frequency transform, to normalize or unnormalize frequencies
            f0 = FT.f0
            BW = FT.BW
            
            # Frequency sampling points
            NpassBand = int(Nsamples/2)
            Noutband   = int((Nsamples+1)/2)
            delta = BW/NpassBand  # Passband sampling step
            freqInband = np.linspace(f0-BW/2+delta/2, f0+BW/2-delta/2, NpassBand)   # Uniform sampling with 1st and last samples at 1/2 step (delta/2) from band transition
            freqOutBandLower = np.linspace(minFreq, f0-BW/2, int(Noutband/2) )
            freqOutBandUpper = np.linspace(f0+BW/2, maxFreq, int((Noutband+1)/2) )
            freq = np.concatenate([freqOutBandLower, freqInband, freqOutBandUpper])
            sEval = 1j * FT.normFreq(freq)
            
            return Sparameters(P, CP, FT, freq)
        # END: Helper function
        #!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        
        x0, NonZeroReal, NonZeroImag, numReal, bounds = self.setupOptimTopology(topologyReal, topologyImag, topologySignReal, topologySignImag)
        
        SPref = setupOptimFreq(P, CP, self.FT, P.minFreq, P.maxFreq, Nsamples)
        
        # Minimum values to test
        S11min = np.max( 20*np.log10(np.abs(SPref.S11)) )   - S11DR  # Minimum value of S11 and S22 to test
        S21min   = np.max( 20*np.log10(np.abs(SPref.S21)) ) - S21DR    # Minimum value of S21 to test
        
        # Optimization parameters
        args = (SPref, S11min, S21min, NonZeroReal, NonZeroImag, numReal, weights, Q) # Function arguments
                
        print "\nCOUPLING MATRIX OPTIMIZATION:"
        if algor == 'L-BFG-S': NonZeroValues = scop.fmin_l_bfgs_b(self.uniformQOptimize, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter, iprint=1 )[0]
        elif algor == 'TNC':     NonZeroValues = scop.fmin_tnc(self.uniformQOptimize, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter )[0]
        elif algor == 'SLSQP':  NonZeroValues = scop.fmin_slsqp(self.uniformQOptimize, x0, fprime=None, args=args,  bounds=bounds, iprint = 2, iter=maxIter)
        else: raise synthError, "Non valid optimization algorithm %s" % algor
        print ""
        
        MatQOut = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
        self.M = MatQOut.M.copy()   # Update self matrix with the optimized version
        self.Qresonators()
        
        
    ##
    # 
    #         Optimization objective function for uniformQ function. The output matrix is forced symmetric, 
    #         with positive imaginary part out of the diagonal and with uniform Q in all the resonators. 
    #         The imaginary part of the diagonal elements is determined so that the Q is constant,so it is not a variable to optimize.
    #         
    #         @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero. 
    #                                             It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
    #         @param SPref = SParameters class instance containing [S] computed from the characteristic polynomials. It is used as the reference to compare.
    #         @param S11min = Minimum value of S11 and S22 to test.
    #         @param S21min  = Minimum value of S21 outside to test.
    #         @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
    #         @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix. 
    #                              The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
    #         @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
    #         @param weights = Optimization weights for the difference in S11(dB), S21(dB) and GD(inband) between the computation from the coupling matrix and from polynomials: 
    #                          List [ weightS11,  weightS22, weightS21, weightGD ]
    #         @param Q =  Prescribed uniform Q for the resonators.
    #         @return goal = Objective function value to minimize.
    #         
    def uniformQOptimize(self, NonZeroValues, SPref, S11min, S21min, NonZeroReal, NonZeroImag, numReal, weights, Q):
        weightS11 = weights[0]
        weightS22 = weights[1]
        weightS21 = weights[2]
        weightGD = weights[3]
        
        # Create an empty matrix and fill it with the non-zero values
        MatQtest = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
        
        # Computation of [S] from the testing coupling matrix
        SPref.fromCouplingMatrix(MatQtest,  checkPassive=False)
        S11      = 20*np.log10(np.abs(SPref.S11))
        S11M   = 20*np.log10(np.abs(SPref.S11M))
        S22      = 20*np.log10(np.abs(SPref.S22))
        S22M   = 20*np.log10(np.abs(SPref.S22M))
        S21   = 20*np.log10(np.abs(SPref.S21))
        S21M = 20*np.log10(np.abs(SPref.S21M))
        
        # Set to minimum values lower than minimum
        S11[S11 < S11min] = S11min
        S11M[S11M < S11min] = S11min
        S22[S22 < S11min] = S11min
        S22M[S22M < S11min] = S11min
        S21[S21 < S21min] = S21min
        S21M[S21M < S21min] = S21min
        
        # Error between [S] from Cp and [S] from CM
        errS11 = np.abs(S11 - S11M)
        errS22 = np.abs(S22 - S22M)
        errS21 = np.abs( S21 -  S21M)
        errGD = np.abs( SPref.groupDelay2 - SPref.groupDelayM ) / np.mean(np.abs(SPref.groupDelay2))    # Normalize to remove 1e-9 scale factor
        
        # Objective function equaling S11 with S11M and S21 with S21M
        goalS11    = np.mean(errS11)
        goalS22    = np.mean(errS22)
        goalS21 = np.mean(errS21)
        goalGD = np.mean(errGD)
        goal = weightS11*goalS11 + weightS22*goalS22 + weightS21*goalS21 + weightGD*goalGD
        
        return goal
        
        
    ##
    # 
    #         Optimize the coupling matrix to achieve a prescribed uniform Q and [S] parameters that comply with a specification mask.
    #         The non-zero coupling matrix entries to optimize are specified by topologyReal and topologyImag parameters, 
    #         which are set by the matrix topology widget in the GUI.
    #         Only the magnitude of [S] is optimized to comply with the S11 and S21 masks.
    #         Group delay in the mask file, if any, is not used in optimization, 
    #         
    #         @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #         @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
    #         @param weights = Optimization weights for the difference in S11 and S21 between the computation from the coupling matrix and from polynomials: 
    #                          List [ weightS11,  weightS22, weightS21, weightGD ]. The last element, weightGD, is not used.
    #         @param Q = Prescribed uniform Q for the resonators.
    #         @param algor = Constrained optimization algorithm. String equal to 'L-BFG-S', 'TNC' or 'SLSQP'. See algorithms description MatrixQ.uniformQ method documentation.
    #         @param maxIter = Maximum number of iterations, or function evaluations, depending on algorithm.
    #         @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between P.minFreq and P.maxFreq.
    #         @param SM = SpecMask class instance, with the mask to compare with.
    #         @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
    #         @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part. 
    #         @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
    #         @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
    #         
    def uniformQMask(self, P, CP, weights, Q, algor, maxIter, Nsamples, SM, topologyReal, topologyImag, topologySignReal, topologySignImag):
        
        x0, NonZeroReal, NonZeroImag, numReal, bounds = self.setupOptimTopology(topologyReal, topologyImag, topologySignReal, topologySignImag)

        # Sparameters class instance, prepared to store [S] computed from CM.
        SPref = Sparameters(P, CP, self.FT)
        
        # Mask [S] parameters to compare with
        SM.specToMask(SPref)    # Needs SP, to get reference values of max(S21) and mean(groupDelay) for relative mask variations. The reference values are obtained from the CP.
        SM.maskSamplesOptimization(self.FT, Nsamples)
        
        # Optimization parameters
        args = (SPref, SM, NonZeroReal, NonZeroImag, numReal, weights, Q) # Function arguments
                
        print "\nCOUPLING MATRIX OPTIMIZATION:"
        if algor == 'L-BFG-S': NonZeroValues = scop.fmin_l_bfgs_b(self.uniformQOptimizeMask, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter, iprint=1 )[0]
        elif algor == 'TNC':     NonZeroValues = scop.fmin_tnc(self.uniformQOptimizeMask, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter )[0]
        elif algor == 'SLSQP':  NonZeroValues = scop.fmin_slsqp(self.uniformQOptimizeMask, x0, fprime=None, args=args,  bounds=bounds, iprint = 2, iter=maxIter)
        else: raise synthError, "Non valid optimization algorithm %s" % algor
        print ""
        
        MatQOut = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
        self.M = MatQOut.M.copy()   # Update self matrix with the optimized version
        self.Qresonators()

        
    ##
    # 
    #         Optimization objective function for uniformQMask function. The output matrix is forced symmetric, 
    #         with positive imaginary part out of the diagonal and with uniform Q in all the resonators. 
    #         The imaginary part of the diagonal elements is determined so that the Q is constant,so it is not a variable to optimize.
    #         Only the magnitude of [S] is optimized to comply with the S11 and S21 masks.
    #         Group delay in the mask file, if any, is not used in optimization, 
    #         
    #         @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero. 
    #                                             It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
    #         @param SPref = SParameters class instance containing method and storage to compute [S] from CM.
    #         @param SM = SpecMask class instance, with the mask to compare with.
    #         @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
    #         @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix. 
    #                              The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
    #         @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
    #         @param weights = Optimization weights for the difference in S11(dB), S21(dB) and GD(inband) between the computation from the coupling matrix and from polynomials: 
    #                          List [ weightS11,  weightS22, weightS21, weightGD ]. The last element, weightGD, is not used.
    #         @param Q = Prescribed uniform Q for the resonators.
    #         @return goal = Objective function value to minimize.
    #         
    def uniformQOptimizeMask(self, NonZeroValues, SPref, SM, NonZeroReal, NonZeroImag, numReal, weights, Q):
        
        weightS11 = weights[0]
        weightS22 = weights[1]
        weightS21 = weights[2]
        
        # Create an empty matrix and fill it with the non-zero values
        MatQtest = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
        
        # Computation of [S] from the testing coupling matrix
        SPref.fromCouplingMatrix(MatQtest, evalS = SM.evalSoptim, checkPassive=False)
        S11M              = 20*np.log10(np.abs( SPref.S11M[:SM.NfreqInBandOptim] ))
        S22M              = 20*np.log10(np.abs( SPref.S22M[:SM.NfreqInBandOptim] ))
        S21MinBand   = 20*np.log10(np.abs( SPref.S21M[:SM.NfreqInBandOptim] ))
        S21MoutBand = 20*np.log10(np.abs( SPref.S21M[SM.NfreqInBandOptim:] ))
        
        # Error between [S] from Mask and [S] from CM
        errS11 = S11M - SM.S11optim
        errS11[errS11<0] = 0
        
        errS22 = S22M - SM.S11optim
        errS22[errS22<0] = 0
        
        errS21MinBand   =  SM.S21inBandOptim - S21MinBand
        errS21MinBand[errS21MinBand<0] = 0
        
        errS21MoutBand = S21MoutBand - SM.S21outBandOptim
        errS21MoutBand[errS21MoutBand<0] = 0
        
        errS21 = np.concatenate([ errS21MinBand, errS21MoutBand ])

        # Objective function
        goalS11    = np.mean(errS11)
        goalS22    = np.mean(errS22)
        goalS21 = np.mean(errS21)
        goal = weightS11*goalS11 + weightS22*goalS22 + weightS21*goalS21
        
        return goal
        
            
    ##
    # 
    #         Save this MatrixQ instance to an ascii file, including energy and power.
    #         Uses self.__str__() to generate the ascii representation of the numpy matrix self.M.
    #         
    #         @param fileName = Name of file, including extension (string).
    #         @param P = Parameter class instance.
    #         @param SP = SParameters class instance.
    #         @param precCM   = Number of significant digits to save in coupling matrix and Q (int).
    #         @param precEP   = Number of significant digits to save in energy and power (int).
    #         
    def saveMat(self, fileName,  P,  SP, precCM, precEP):
        #fout = codecs.open(fileName, 'w',  'utf-8')    # codecs.open does not translate /n into /r/n in windows.
        fout = open(fileName, 'w')
        fout.write(u'# %s; version %s\n\n' % pkgNameVersion(__doc__) )
        fout.write(u'# Coupling matrix name\n')
        fout.write(u'Name: %s\n\n' % self.name)
        fout.write(u'# Number of non-resonant nodes at the source and load sides. Includes the source and load nodes.\n')
        fout.write(u'extraNodesS: %d\n' % self.extraNodesS)
        fout.write(u'extraNodesL: %d\n\n' % self.extraNodesL)
        
        fout.write(u'# Vector with node scaling factors\n')
        st = u"nodeScalingFactor: %s\n\n" %  ["%.4g" % elem for elem in self.nodeScalingFactor]
        fout.write(st.replace("'", ""))
        
        fout.write(u'# Flag which indicates if it is allowed to do a matrix rotation (impossible after a resonant node scaling)\n')
        fout.write(u'flagRotateAllowed: %s\n\n' % self.flagRotateAllowed)
        
        fout.write(u'# Quality factor of resonators\n')
        st = u"Q: %s\n\n" %  ["%.4g" % elem if abs(elem) != np.inf else "%s" % str(elem) for elem in self.Q]
        fout.write(st.replace("'", ""))
        
        fout.write(u'# Coupling matrix (each line is a matrix row):\n')
        st = self.__str__(precCM)
        fout.write(st)
        
        # Compute energy and power for this matrix
        SP.energyPowerFromCM(self)
        
        # Row and column labels
        size = self.M.shape[0]
        first = self.extraNodesS
        last = size - self.extraNodesL
        labels = ['S']
        for n in range(1, first ): labels.append('NR' + str(n))
        labels.extend(  [ 'R'+str(elem) for elem in range(1,  last-first+1) ] )
        for n in range(first,  first+self.extraNodesL-1 ): labels.append('NR' + str(n))
        labels.append('L')
        
        powerRes = SP.powerRes
        energyRes = SP.energyRes
        powerResTotal = powerRes.sum(axis=1)
        energyResTotal = energyRes.sum(axis=1)
        powerCoup =SP.powerCoup
        energyCoup = SP.energyCoup
        
        # Set column width
        colWidth = 10 + precEP
        
        Nf = np.size(energyRes,0)    # Number of frequencies
        Nr = np.size(energyRes,1)    # Number of resonators
        
        fout.write ('\n# ENERGY AND POWER\n')
        
        fout.write ('\n# Stored energy in the resonators:\n')
        st = '%-*s' % (colWidth, '#GHZ')
        for n in range(Nr): st += 'R%-*d' % (colWidth-1, n+1)
        st += 'Total\n'
        fout.write(st)        
        st = ''
        for m in range(Nf):
            freq = str(SP.freq[m] / 1e9)
            st += '%-*s' % (colWidth, freq)
            for n in range(Nr): 
                stNum = '%.*g' % (precEP, energyRes[m, n])
                st += '%-*s' % (colWidth,  stNum) 
            stNum = '%.*g' % (precEP, energyResTotal[m])
            st += '%-*s\n' % (colWidth,  stNum)
        fout.write(st) 
        
        fout.write ('\n# Dissipated power in the resonators:\n')
        st = '%-*s' % (colWidth, '#GHZ')
        for n in range(Nr): st += 'R%-*d' % (colWidth-1, n+1)
        st += 'Total\n'
        fout.write(st)        
        st = ''
        for m in range(Nf):
            freq = str(SP.freq[m] / 1e9)
            st += '%-*s' % (colWidth, freq)
            for n in range(Nr): 
                stNum = '%.*g' % (precEP, powerRes[m, n])
                if stNum == '-0': stNum = '0'
                st += '%-*s' % (colWidth,  stNum) 
            stNum = '%.*g' % (precEP, powerResTotal[m])
            st += '%-*s\n' % (colWidth,  stNum)
        fout.write(st)
        
        if len(powerCoup) > 0:
            fout.write ('\n# Dissipated power in the couplings:\n')
            st = '%-*s' % (colWidth, '#GHZ')
            for coup in powerCoup: 
                (m, n) = coup[0]
                elem = '%s-%s' % (labels[m], labels[n]) 
                st +='%-*s' % ( colWidth, elem)
            st += '\n'
            fout.write(st)        
            st = ''
            for m in range(Nf):
                freq = str(SP.freq[m] / 1e9)
                st += '%-*s' % (colWidth, freq)
                for coup in powerCoup:
                    stNum = '%.*g' % (precEP, coup[2][m])
                    if stNum == '-0': stNum = '0'
                    st += '%-*s' % (colWidth,  stNum) 
                st += '\n'
            fout.write(st)
        
        fout.close()
        
        
    ##
    # 
    #         Convert numpy matrix to string representation.
    #         Rounds coupling matrix elements to given precision using complexStr() function.
    #         @param prec     = Number of significant digits (int). Default 5.
    #         @return st = String representation of the numpy matrix self.M.
    #         
    def __str__(self,  prec = 5):

        M = self.M
        (mrows, ncols) = M.shape
        
        # Set column width
        colWidth = 5                                                                           # Blank space between columns
        if (abs(M.real)>10**-prec).any():  colWidth += prec+4       # Space for the real part
        if (abs(M.imag)>10**-prec).any(): colWidth += prec+5       # Space for the imaginary part
        
        st = '# ' + self.name + u': Low-pass matrix (normalised)\n'
        for m in range(mrows):
            for n in range(ncols):
                tmp = complexStr(M[m, n],  prec)
                st += u'%-*s' % (colWidth,  tmp)
            st += u'\n'
            
        self.lp2bp2cbw()
        
        st += '\n\n# ' + self.name + u': Coupling bandwidth matrix (MHz)\n'
        for m in range(mrows):
            for n in range(ncols):
                tmp = complexStr(self.Mcbw[m, n],  prec)
                st += u'%-*s' % (colWidth,  tmp)
            st += u'\n'
        
        st += '\n\n# ' + self.name + u': Coupling bandwidth matrix with transmission lines (MHz)\n'
        for m in range(mrows):
            for n in range(ncols):
                tmp = complexStr(self.McbwTL[m, n],  prec)
                st += u'%-*s' % (colWidth,  tmp)
            st += u'\n'
            
        return (st)
        
        
# End class: Sparameters
#!!!!!!!!!!!!!!!!!!!!!!!!!#


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin class: CouplingMatrices

##
# 
#     Coupling matrices class.
#     If \f$ [Z] \f$ is the impedance matriz that relates the current at the network loops with the voltage sources within these loops,
#     the compupling matriz \f$ M \f$ is such that [TN 102.1, sec. 3.2]:
#     \f[
#         [Z] = [jM + sI]
#     \f]
#     @todo Call numpy functions as methods of array class, whenever possible.
#     
class CouplingMatrices(object):
    
    ##
    # 
    #         Compute Coupling Matrices for all the selected topologies.
    #         @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #         @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
    #         @param Yo = Termination admitance.
    #         @param FT = Frequency transformation class instance.
    #         @param SP = SParameters class instance.
    #         @return CM = CouplingMatrices class instance. 
    #         
    def __init__(self, P,  CP,  Yo,  FT, SP):
        
        ##
        # Frequency transform used to compute these coupling matrices. 
        # It must be used when creating new matrices or reading matrices from disk.
        self.FT = FT
        
        ##
        # SParameters class instance containing the frequency sampling and SP.fromCouplingMatrix() method.
        # Will be used in the sensitivity analysis.
        self.SP = SP
        
        ##
        # Uniform Q folded coupling matrix N+4.
        self.fCM_uniQ = None
        
        ##
        # List of available coupling matrices
        self.listM = [ ]
        
        ##
        # Last index of CM selection comboBox. Initially it is equal to zero(first matrix in list).
        self.indexCM = 0
        
        ##
        # Roundoff error in the computation of the Tranversal Coupling Matrix N+2.
        # The error is obtained as the maximum of the largest of the remainders in the computation of \f$ y_{22n}(s) \f$ and  \f$ y_d(s) \f$ 
        # and the constant term in the partial fraction expansion of \f$ y_{21n}(s) / Y_d(s) \f$
        self.CM_error = None
                
        # Pre-computed coupling matrices
        tCM = self.transCouplingMatrix(P,  CP,  Yo, FT)   # Transversal Coupling Matrix N+2. Definition in TN 102.1, sec. 3.2, Fig. 2b), eq. (2), (3) and (4).
        fCM = self.tcm2fcm(tCM)                         # Folded coupling matrix N+2.
        if CP.flagCase1: fCM_uniQ = self.uniformQFCMOrder4(fCM)               # Folded coupling matrix N+4 with uniform Q.
        if CP.flagCase2: fCM_uniQ = self.uniformQFCMOrder6(fCM)               # Folded coupling matrix N+4 with uniform Q.
        if CP.flagCase3: fCM_uniQ = self.uniformQFCMOrder6Case3(fCM)     # Folded coupling matrix N+4 with uniform Q.
        
        # Setting to zero small elements before matrix rotations can be worse because numpy sets a 0+0j where you have set a zero, and NaN appear in trigonometric operations
        #fCM.M [np.abs(fCM.M) < 2*self.CM_error] = 0. # Set to zero small elements of the CM, to avoid problems in matrix rotations

        arrow = self.tcm2arrow(tCM)   # arrow is never equal to None because it can always be computed.
        #arrow.M [np.abs(arrow.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        transComplexZeros = np.array(P.transComplexZeros)
        transComplexZerosWithSym = np.concatenate((transComplexZeros, -transComplexZeros.conj()))
        normImagZeros      = [ 1j*elem for elem in FT.normFreq(np.array(P.transImagZeros)) ]
        normComplexZeros = transComplexZerosWithSym.real + 1j * FT.normFreq(transComplexZerosWithSym.imag)
        TZ = np.concatenate(( normImagZeros, normComplexZeros ))
        if P.genChebyLP: TZ = np.concatenate([TZ, np.array(P.LPcomplexZeros).real + 1j * FT.normFreq(np.array(P.LPcomplexZeros).imag), -np.conjugate(np.array(P.LPcomplexZeros).real + 1j * FT.normFreq(np.array(P.LPcomplexZeros).imag)), np.array(P.LPrealZeros),  -np.array(P.LPrealZeros) ])

        culdesac = self.fcm2culdesac(fCM, TZ) # May be equal to None if it is not posible to compute this matrix
        #if culdesac:
        #    culdesac.M [np.abs(culdesac.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        triplet = self.arrow2triplet(arrow, TZ)  # May be equal to None if it is not posible to compute this matrix
        #if triplet:
        #    triplet.M [np.abs(triplet.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        cascadedQuartets = self.fcm2cqs(fCM, 0, 'None')  # May be equal to None if it is not posible to compute this matrix
        #if cascadedQuartets:
        #    cascadedQuartets.M [np.abs(cascadedQuartets.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        pfitzenmaier = self.fcm2pfitzenmaier(fCM)  # May be equal to None if it is not posible to compute this matrix
        #if pfitzenmaier:
        #    pfitzenmaier.M [np.abs(pfitzenmaier.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        inlineAsym = self.fcm2inlineAsymmetric(fCM)  # May be equal to None if it is not posible to compute this matrix
        #if inlineAsym:
        #    inlineAsym.M [np.abs(inlineAsym.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations

        inlineSym = self.fcm2inlineSymmetric(fCM, TZ)  # May be equal to None if it is not posible to compute this matrix
        #if inlineSym:
        #    inlineSym.M [np.abs(inlineSym.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
        
        # Append pre-computed coupling matrices to the list:
        self.listM.append( tCM)
        self.listM.append( fCM ) 
        if CP.flagCase1 or CP.flagCase2 or CP.flagCase3:
            self.listM.append( fCM_uniQ )
        self.listM.append( arrow )
        if culdesac:           self.listM.append( culdesac )
        if triplet:            self.listM.append( triplet )
        if cascadedQuartets:   self.listM.append( cascadedQuartets )
        if pfitzenmaier:       self.listM.append( pfitzenmaier )
        if inlineAsym:         self.listM.append( inlineAsym )
        if inlineSym:          self.listM.append( inlineSym )
        
        ##
        # Current coupling matrix, after rotations. Is the first in the list.
        self.MatQ = self.listM[0]
   
        ##
        # Backup list of self.MatQ  values, for undo. Must use deep copy, since self.MatQ will be modified by edit actions.
        self.bkpMatQ = [ self.MatQ.copy() ] 
        
        ##
        # Backup list of actions
        self.bkpActions = [ ]


    ##
    # 
    #         Convenience function that computes the threshold to determine if coupling matrix elements (or their real or imaginary parts) are null.
    #         @param st = String with relevant info. Default None.
    #         @return eps = Threshold.
    #         
    def matrixElementsEps(self, st = None):
        eps = self.CM_error*10      # Let threshold be rather close to error estimation, since elements will be filtered later in the Coupling Matrix window by using "Precision" widget.
        if eps < 1e-9: eps = 1e-9
        if eps > 1e-5: eps = 1e-5
        
        if verbose and st:
            myPrint('Truncation threshold = %.3e for CM_error = %.3e (%s)' % (eps, self.CM_error, st) )
        
        return eps
        
    
    ##
    # 
    #         Convenience function that updates the error estimation in self.CM_error
    #         @param CM_error = New error estimation.
    #         @param st = String with relevant info. Default None.
    #         
    def updateCM_error(self, CM_error, st = None):
        
        if verbose and st:
            myPrint('Updating CM error estimation %.2e with %.2e (%s)\n' % (self.CM_error, CM_error, st))
            
        self.CM_error = max([self.CM_error,  CM_error])
            
    
    ##
    # 
    #             Computation of the \f$ N+2 \f$ transversal coupling matrix of the filter defined by its caracteristic polynomials and constants.
    #             Polynomials are of type numpy.poly1d.
    #             This function sets the self.CM_error attribute to the roundoff error in the computation of M, 
    #             given by the largest of the remainders in the computation of \f$ y_{22n}(s) \f$ and  \f$ y_d(s) \f$ 
    #             and the constant term in the partial fraction expansion of \f$ y_{21n}(s) / Y_d(s) \f$
    #             @param P  = Parameters class instance, containing filter, synthesis and sweep parameters.
    #             @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
    #             @param Yo = Termination admitance.
    #             @param FT = Frequency transformation class instance.
    #             @return MatQ = MatrixQ class instance, containing the N+2 transversal coupling matrix.
    #             
    #         \section Smat Scattering matrix [S]
    #         
    #         The first step to find the \f$ N+2 \f$ transversal coupling matrix of a two-port network is finding the associated scattering matrix \f$ [S] \f$. 
    #         
    #         By definition, the scattering matrix of a two-port network is  [TN 102.1, eq (1)]:
    #         \f[
    #         [S]=\begin{bmatrix} S_{11}(s) & S_{12}(s) \\ S_{21}(s) & S_{22}(s) \end{bmatrix} = \frac{1}{S_{d}(s)} \begin{bmatrix} S_{11n}(s) & S_{12n}(s) \\ S_{21n}(s) & S_{22n}(s) \end{bmatrix}
    #         \f]
    #         
    #        The [S] parameters are computed from the polinomials as:
    #        \f[ 
    #             S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)} 
    #         \f]
    #         \f[
    #             S_{21} = S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)}
    #         \f]
    #         \f[
    #             S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)}
    #         \f]
    #         except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where [S] parameters are computed from the lossless case above as:
    #         \f[
    #             S_{11} = k_{11} S_{11\mathrm{lossless}}
    #         \f]
    #         \f[
    #             S_{21} = S_{12} = k_{21} S_{21\mathrm{lossless}}
    #         \f]
    #         \f[
    #             S_{22} = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right) + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11  \mathrm{lossless}}
    #         \f]   
    #         
    #         and hence the numerators and denominator of the \f$ S_{ij} \f$ rational polynomials are:
    #         \f[
    #             S_{11n}(s) = \frac{F(s)}{\epsilon_{R}}
    #         \f]
    #         \f[
    #             S_{12n}(s) = S_{21n} = \frac{P(s)}{\epsilon}
    #         \f]
    #         \f[
    #             S_{d}(s) = E(s)
    #         \f]
    #         \f[
    #             S_{22n}(s) = \frac{F_{22}(s)}{\epsilon_R}
    #         \f]
    #         except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where they are computed as
    #         \f[
    #             S_{11n}(s) = k_{11} \frac{F(s)}{\epsilon_{R}}
    #         \f]
    #         \f[
    #             S_{12n}(s) = S_{21n} = k_{21} \frac{P(s)}{\epsilon}
    #         \f]
    #         \f[
    #             S_{d}(s) = E(s)
    #         \f]
    #         \f[
    #             S_{22n}(s) = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right) \, S_d(s) + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11n}(s)
    #         \f]  
    #         In practice, the software divides \f$ \epsilon_R \f$ by \f$ k_{11} \f$ and \f$ \epsilon \f$ by \f$ k_{21} \f$ , 
    #         and stores the updated values to automatically apply the \f$ k_{11} \f$ and \f$ k_{21} \f$ whenever computing \f$  S_{11n} \f$ and \f$  S_{21n} \f$ .
    #                 
    #         \section Ymat Admitance matrix [Y]
    #     
    #         The next step is the computation of \f$ y_{21}(s)\f$ and \f$ y_{22}(s) \f$ of the admitance matrix \f$ [Y] \f$. 
    #         [Y] is related to [S] as [TN 102.1, eq (12)]:
    #         \f[
    #             [Y]= \frac{1}{y_{d}(s)} \begin{bmatrix} y_{11n}(s) & y_{12n}(s) \\ y_{21n}(s) & y_{22n}(s) \end{bmatrix}
    #         \f]
    #         Where:
    #         \f[
    #             y_{21n}(s)= -2 \, S_{21n}(s)
    #         \f]
    #         \f[
    #             y_{22n}(s)=\frac{(S_{d}(s)+S_{11n}(s))(S_{d}(s)-S_{22n}(s))+S_{12n}S_{21n}}{S_{d}(s)}
    #         \f]
    #         \f[
    #             y_{11n}(s)=\frac{(S_{d}(s)-S_{11n}(s))(S_{d}(s)+S_{22n}(s))+S_{12n}S_{21n}}{S_{d}(s)}
    #         \f]
    #         \f[
    #             y_{d}(s)=\frac{(S_{d}(s)+S_{11n}(s))(S_{d}(s)+S_{22n}(s))-S_{12n}S_{21n}}{Y_{0}(s)S_{d}(s)}
    #         \f]
    #         \f$ Y_{0} \f$ = Load admitance at the end of the network.
    #         
    #         \section Mmat N+2 Coupling Matrix
    #         A partial fraction expansion of \f$ [y] \f$ yields [Cameron, eq. (8.41)] [TN 101.1, eq (12)]:
    #         \f[
    #         [Y]= \frac{1}{y_{d}(s)} \begin{bmatrix} y_{11n}(s) & y_{12n}(s) \\ y_{21n}(s) & y_{22n}(s) \end{bmatrix}
    #         = j \begin{bmatrix} 0 & K_{\infty} \\ K_{\infty} & 0 \end{bmatrix} 
    #         + \sum_{k=1}^N \frac{1}{s-j\lambda_k} \cdot  \begin{bmatrix} r_{11k}(s) & r_{12k}(s) \\ r_{21k}(s) & r_{22k}(s) \end{bmatrix}
    #         + \begin{bmatrix} KK_{11} & 0 \\ 0 & KK_{22} \end{bmatrix} 
    #         \f]
    #      
    #         <ul>
    #         
    #         <li> Computation of the residuals \f$ r_{22k} \f$ and \f$ r_{21k} \f$:
    #      
    #         The residuals \f$ r_{22k} \f$ and \f$ r_{21k} \f$ can be found as [Cameron, pp. 294, footnote]:
    #         \f[
    #             r_{22k}(s) = \left.\frac{y_{22n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}} \qquad
    #             r_{11k}(s) = \left.\frac{y_{11n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}} \qquad
    #             r_{12k}(s) = r_{21k}(s) = \left.\frac{y_{21n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}}
    #         \f]
    #         
    #         where:
    #             <ul>
    #             <li> \f$ y'_{d}(s) \f$ is the first derivative of \f$ y_{d}(s)\f$. It is variable Ydp in the code.
    #             <li> \f$ j\lambda_{k} \f$ are the roots of  \f$ y'_{d}(s) \f$.  
    #             </ul>
    #         The residuals are actually computed using scipy.signal.residue() function in order to obtain the constant terms:
    #             <ul>
    #             <li>In \f$ [Y_{21}] \f$ the constant term is the source-load coupling term \f$ j K_{\infty} \f$ :
    #             \f[ \frac{y_{21n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{21k}}{s-j\lambda_k} + j K_{\infty}
    #             \f]
    #             <li>In \f$ [Y_{22}] \f$ and \f$ [Y_{11}] \f$ , "Prescribed Insertion Loss technique" case 1, there is a constant term \f$ \text{\it{KK}} \not= 0 \f$ 
    #             that must be included in the N+2 transversal coupling matrix  [TN 102.1, eq (13)]:
    #             \f[ \frac{y_{22n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{22k}}{s-j\lambda_k} + KK_{22}
    #             \f]
    #             
    #             \f[ \frac{y_{11n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{11k}}{s-j\lambda_k} + KK_{11}
    #             \f]
    #             </ul>
    #         
    #         <li> Computation of \f$ K_{\infty} \f$:
    #      
    #         \f$ K_{\infty} \f$ is a real constant equal to 0 except for fully canonical filters 
    #         where the number of finite-position transmission zeros  \f$ n_{fz} \f$ is equal to the filter degree \f$ N \f$. 
    #         For lossless filters it can be computed as [Cameron, eq (8.43a)]:
    #         \f[
    #             K_{\infty} = \frac{\epsilon_{R}}{\epsilon} \frac{1}{(\epsilon_{R}+1)}
    #         \f]
    #         However, this formula is not valid for lossy synthesis, so we set \f$ K_{\infty} \f$ equal to the imaginary part of the constant term in the polynomial expansion of \f$ y_{21n}'(s) / y_d(s) \f$ .
    #         
    #         <li> Computation of the eigenvectors \f$ T_{1k} \f$ and \f$ T_{Nk} \f$:
    #          
    #         The eigenvectors \f$ T_{Nk} \f$ and \f$ T_{1k} \f$ can be found as [Cameron eq. (8.55)]:
    #         \f[
    #             T_{Nk} = \sqrt{r_{22k}}
    #         \f]
    #         \f[
    #             T_{1k} = r_{21k} / \sqrt{r_{22k}}
    #         \f]
    #         
    #         <li> Build the \f$ N+2 \f$ transversal coupling matrix \f$ [M] \f$:
    #      
    #         The \f$ N+2 \f$ transversal coupling matrix is built from the orthogonal vectors \f$ T_{Nk} \f$ and \f$ T_{1k} \f$ 
    #         and from the eigenvalues \f$ \lambda_{k} \f$, as follows [Cameron, fig 8.18]:
    #         \f[
    #             [M] = \begin{bmatrix}
    #                     -j KK_{11}       & \cdots    & T_{1k}        & \cdots    & K_{\infty} \\ 
    #                     \vdots      & \ddots    &               & 0         & \vdots    \\
    #                     T_{1k}      &           & -\lambda_{k}  &           & T_{Nk}    \\
    #                     \vdots      & 0         &               & \ddots    & \vdots    \\
    #                     K_{\infty}  & \cdots    & T_{Nk}        & \cdots    & -j KK_{22} 
    #                 \end{bmatrix}
    #         \f]
    #         The steps to build the matrix are the following:
    #             <ol>
    #             <li> We put the vector T1k at the first row (between the first and the last column not-included)
    #             <li> We put the value of Kinf at the last column of the first row
    #             <li> We put the vector TNk at the last column (between the first and the last row not-included)
    #             <li> As the matrix M is symmetric, we copy the upper triangle to the lower
    #             <li> We put the lambdak (with opposite sign) at the main diagonal 
    #                 (between the first and the last element not-included)
    #             <li> If there is a remaining constants \f$ KK_{11} \f$ and  \f$ KK_{22} \f$
    #                  in the partial fraction expansion of \f$ Y_{11n}/Y_d \f$ and \f$ Y_{22n}/Y_d \f$, 
    #                  set M[0,0] to \f$ -j KK_{11} \f$ and M[N+1,N+1] to \f$ -j KK_{22} \f$ .
    #             </ol>
    #         </ul>
    #         @todo The function "r21k,  p21k,  remainder = sc.signal.residue(Y21np,Yd)" does not give the correct remainder when it is complex but in the tested cases it does not afect. This function should be improved for the future.
    #         
    def transCouplingMatrix(self, P,  CP,  Yo, FT):
        
        # Degree of the filtering function
        N = CP.Es.o
        # Number of finite-position transmission zeros (TZs) that have been prescribed.
        nfz = CP.Ps.o
        
        # Scattering matrix [S] [TN 102.1, eq (1)]
        S11n = CP.Fs/CP.epsR
        S12n = CP.Ps/CP.eps
        S21n = CP.Ps/CP.eps
        Sd = CP.Es
        
        if P.synthTech.title() == u'Prescribed Insertion Loss' and P.nuqTech == u'kS21+kS11':
            k11 = 10**(-P.nuqK11c1/20.0)  # Insertion lossess
            k22 = 10**(-P.nuqK22c1/20.0)  # Insertion lossess
            k21 = 10**(-P.nuqK21c1/20.0)  # Insertion lossess
            S22n = k22 * ( 1- (k21/k11)**2) * Sd + (k21/k11)**2 * S11n
        else:
            S22n = CP.F22s/CP.epsR
       
        # Admitance matrix [Y] [TN 102.1, eq (12)]:
        Y21n = -2*S21n
        
        #  Computation of Y22n and Yd with polynomial division. 
        # This function does not work well with the example in "Prescribed_insertion_losses_K11+k21_asymmetrical.par"
        # (the example in TN501), while numpy.polydiv from numpy v1.3.0 works perfectly.
#        Y22n,  error1 = self.polyExDiv( (Sd+S11n)*(Sd-S22n) + S12n*S21n,  Sd)
#        Y11n,  error2 = self.polyExDiv( (Sd-S11n)*(Sd+S22n) + S12n*S21n,  Sd)
#        Yd, error3      = self.polyExDiv( (Sd+S11n)*(Sd+S22n) - S12n*S21n,  Sd)
        
        Y22n,  error1 = ((Sd+S11n)*(Sd-S22n) + S12n*S21n)  /  Sd
        Y11n,  error2 = ((Sd-S11n)*(Sd+S22n) + S12n*S21n)  /  Sd
        Yd, error3      = ((Sd+S11n)*(Sd+S22n) - S12n*S21n)  /  Sd
        # Although the round-off error is larger than with polynomial evaluation and polyfit, it is given by the division remainder.
        # numpy.polydiv and poly1d division operator fail (for numpy 1.1) with complex polynomial coeeficients.
        
        if Yd.order != N: raise synthError,  'Transversal Coupling Matrix computation: Order of Yd(s) different than order of E(s).'
        CM_error = abs(np.concatenate( (error1.coef, error2.coef, error3.coef),  1) ).max()
        
        if verbose: myPrint('\nTCM (N+2): error in polinomial division (Y22n, Y11n, Yd) = (%.2e, %.2e, %.2e)' % (abs(error1.coef).max(), abs(error2.coef).max(), abs(error3.coef).max() ))
        
        # Normalization to have highest order coefficients equal to 1
        Y21n = Y21n/Yd.coef[0]
        Y22n = Y22n/Yd.coef[0]
        Y11n = Y11n/Yd.coef[0]
        Yd = Yd/Yd.coef[0]
    
        # Computation of the residuals r22k and r21k  [TN 102.1, eq (13)], [Cameron, pp. 294]
        r22k,  p22k, KK22  = sc.signal.residue(Y22n,Yd)
        r11k,  p11k, KK11  = sc.signal.residue(Y11n,Yd)
        #r21k,  p21k, KK2 = sc.signal.residue(Y21np,Yd)
        r21k,  p21k, KK21 = sc.signal.residue(Y21n,Yd)
        lambdak = p21k / 1j        
        
        # Kinf [Cameron, eq (8.43a)]
        if N<>nfz: # For non canonical filters
            Kinf = 0
        else: # For canonical filters
            #Kinf = CP.epsR/(CP.eps*(CP.epsR+1.)) # Only for lossless filters, does not work with lossy synthesis
            Kinf = KK21[-1].imag # KK22 is a polynomial (degree 0 if there are no round-off errors) and KK22[-1] is the constant term
            
        # Recomputation of Y21n for canonical filters [Cameron, eq. (8.44)]. If the filter is not canonical Y21n = Y21np
        # It is not necessary, since we already have the correct residues.
        #Y21np = Y21n - 1j*Kinf*Yd
        #r21k,  p21k, KK2 = sc.signal.residue(Y21np,Yd)
        
        # Computation of the eigenvectors T1k and TNk [Cameron eq. (8.55)]:
        iLambdak = np.real(lambdak).argsort()      # Sorting order: Decreasing real part of lambdak
        lambdak = lambdak[iLambdak]
        TNk = np.sqrt(r22k[iLambdak])
        T1k = r21k[iLambdak]/np.sqrt(r22k[iLambdak])
        
        # Build the N+2 transversal coupling matrix [M] [Cameron, fig 8.18]:
        M = 1j*np.zeros((N+2,N+2))
        
        # Discard constant term smaller than roundoff error
        if abs(KK11) <= CM_error: KK11 = np.array([0])
        if abs(KK22) <= CM_error: KK22 = np.array([0])
        
        M[0,1:N+1] = T1k.T          # 1.- We put the vector T1k at the first row (between the first and the last column not-included).
        M[0,N+1] = Kinf                # 2.- We put the value of Kinf at the last column of the first row. 
        M[1:N+1,N+1] = TNk        # 3.- We put the vector TNk at the last column (between the first and the last row not-included)
        M = M + M.T                      # 4.- As the matrix M is symmetric, we copy the upper side to the lower
        # 5.- We put the lambdak (with opposite sign) at the main diagonal (between the first and the last element not-included)
        # with the remaining constants KK11, KK22 in the partial fraction expansion of Y22n/Yd, in M[0,0] and M[N+1,N+1]
        M = M - np.diag(np.concatenate((1j*KK11, lambdak, 1j*KK22), 1))
        
        # If synthesis is lossless, in theory the imaginary part of couplings is roundoff error.
        # However, there are cases such as "Cameron - Example TCM FCM sec 8.4.3.par" in which 
        # the residuals have imaginary parts of the order 1e-5, which are not numerical errors and are necessary
        # to make [S] parameters computed from CM agree [S] parameters computed from the CP.
        
        if (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+pole/zero'):
            # Remove the column and row corresponding to the (s-a)/(s+a) zero/pole that has been added in case 3 of Prescribed Insertion Loss synthesis.
            assert CP.nuqDelta is not None

            # Old way to remove row and column (Txema)
#            eps = 1e-6
#            ind = np.nonzero( np.abs( - np.imag(np.diag(M)) / CP.nuqDelta - 1 ) > eps )  
#            M = M[np.ix_(ind[0], ind[0])]
            
            # New approach (JMR): Find the minimum in order to avoid problems with thresholding
            # We could get two poles (or no pole) closer than eps to CP.nuqDelta
            # It is better to get the closest pole to CP.nuqDelta
            indDelta = np.argmin( np.abs(np.imag(np.diag(M)) + CP.nuqDelta) )
            M = np.delete(np.delete(M.copy(), indDelta, 1), indDelta, 0)
            
            
        # If Minimum Insertion Loss technique with finite Q, put losses in the diagonal of M
        if P.synthTech.title() == u'Minimum Insertion Loss' and P.finiteQo:
            for n in range(N):
                M[n+1, n+1] += -1j  / (P.Qo * FT.FBW)
                
        if P.synthTech.title() == u'Predistortion':
            for n in range(N):
                M[n+1, n+1] += -1j  / (P.QoPredis * FT.FBW)
            
        # Set error in computation of coupling matrix
        self.CM_error = CM_error
        if verbose: myPrint('TCM (N+2): overall error = %.2e\n' % CM_error)
        
        MatQ = MatrixQ(self, M,  'Transversal (N+2)', 1,  1, FT,  fileExt = '.tcm')
        
        # Make the sign of coupling (S,R1) equal to sign of (L,R1). 
        # This is equivalent to changing sign of poly P(s)
        # This is necessary to obtain the same TCM and FCM for Cameron example at section 8.4.3 and software CMS Demo.
        if np.sign(MatQ.M[1,0]) != np.sign(MatQ.M[1,1]): MatQ.scaleNode(0, -1, True)
        
        return MatQ
     
       
    ##
    # 
    #         Polynomial exact division.
    #         Computes quot and remain polynomials so that num(s) = quot(s)*den(s) + error(s).
    #         @param  num   = Dividend.
    #         @param  den   = Divisor.
    #         @return quot  = Quotient.
    #         @return error = Error in the computation of the quotient. Is equal to: num(s)-quot(s)*den(s)
    #         
    #         \section algPolyDiv Algorithm
    #         Since numpy.polydiv or the equivalent numpy.poly1d polynomial division fails for complex coefficients, 
    #         at least in numpy v1.1.1, here we evaluate num(s)/den(s) for a finite number of points in the imaginary axis
    #         and use polyfit to build the polynomial that passes trhough this points. The remainder is the error.
    #         
    #         This function does not work well with the example in "Prescribed_insertion_losses_K11+k21_asymmetrical.par"
    #         (the example in journal paper publication), while numpy.polydiv from numpy v1.3.0 works perfectly.
    #         
    #         
    def polyExDiv(self,  num,  den):
        N = num.o - den.o       # Order of Quotient
        
        sValues = 1j * np.linspace(-2.0,  2.0,  3*(N+1))        # Use 3x the minimum number of points needed. Same result as using 40 points always.
        quot = np.poly1d(np.polyfit( sValues,  num(sValues)/den(sValues),  num.o-den.o ))
        error = num - quot*den
        return quot,  error
        
        
    ##
    # 
    #         This function transforms the transversal \f$N+2\f$ coupling matrix \f$M\f$ into the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$.
    #         @param MatQ = MatrixQ class instance, containing the transversal \f$N+2\f$ coupling matrix.
    #         @return MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         
    #         [TN 102.2 sec. 9.2]
    #         The procedure to obtain the matrix [FCM] consists in setting to zero step by step elements of the matrix [M] using the following:
    #         
    #         -Step 1: We remove all the elements of the row 1 (In total there are nerk(1) elements).
    #         
    #         -Step 2: We remove all the elements of the column 1 (In total there are neck(1) elements).
    #         
    #         -Step 3: The same as in step 1 with the row 2 (In total there are nerk(2) elements).
    #         
    #         -Step 4: The same as in step 2 with the column 2 (In total there are neck(2) elements).
    #         
    #         \f$\vdots\f$
    #         
    #         -Step numSteps-1:  We remove all the elements of the row numRows (In total there are nerk(numRows) elements).
    #         
    #         -Step numSteps: We remove all the elements of the column numColumns (In total there are nerk(numColumns) elements).
    #         
    #         The last two steps consist in removing the last remaining row and the last remaining column of the matrix (it is not necessary that they are in that order).                                                                                                                                                                
    #         
    def tcm2fcm(self, MatQ):
        Mfcm = MatQ.copy()
        Mfcm.setName('Folded (N+2)')
        Mfcm.setFileExt('.fcm')
        itr=0
        itc=0
        
        # Size of the matrix
        N = np.size(Mfcm.M,0)
        
        # Total number of steps
        numSteps = (N-3)*(N-2)/2
    
        # Total number of rows numRows and columns numColumns to be manipulated
        if N%2:
            numRows = (N-3)/2
            numColumns = (N-3)/2
        else:
            numRows = (N-2)/2
            numColumns = (N-4)/2
        for k in range(numRows+numColumns):
            if (k+1)%2:
                # Elimination process of one row
                nerk = N-3-2*itr
                for l in range(nerk):
                    m = itr                                                                                 
                    n = N-l-itr-2
                    i = N-l-itr-3
                    j = N-l-itr-2
                    rotAng = Mfcm.rotAngleEliminate('trigonometric', i, j, m, n)
                    Mfcm.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
                itr = itr+1
            else:
                # Elimination process of one column
                neck = N-4-2*itc
                for l in range(neck):
                    m = 2+l+itc
                    n = N-itc-1
                    i = 2+l+itc
                    j = 3+l+itc
                    rotAng = Mfcm.rotAngleEliminate('trigonometric', i, j, m, n)
                    Mfcm.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
                itc = itc+1
        Mfcm.Qresonators()
        return Mfcm
        
    ##
    # 
    #         This function transforms the transversal \f$N+2\f$ coupling matrix \f$M\f$ into the coupling matrix in Arrow or Wheel Form \f$Marrow\f$.
    #         @param MatQ = MatrixQ class instance, containing the transversal \f$N+2\f$ coupling matrix.
    #         @return MatQarrow =  MatrixQ class instance, containing the Arrow topology coupling matrix. 
    #         
    #         [Cameron]
    #         
    def tcm2arrow(self, MatQ):
    
        Marrow = MatQ.copy()
        Marrow.setName('Arrow (N+2)')
        Marrow.setFileExt('.arrow')
        
        # Size of the matrix
        N = np.size(Marrow.M,0)-2
        
        for i in range(1, N):
            for j in range(i+1, N+1):
                rotAng = Marrow.rotAngleEliminate('trigonometric', i, j, i-1, j)
                Marrow.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        Marrow.Qresonators()
        return Marrow

    ##
    # 
    #         This function transforms the coupling matrix in Arrow or Wheel Form \f$Marrow\f$ into the coupling matrix in Trisections Form \f$Mtriplet\f$.
    #         @param MatQarrow = MatrixQ class instance, containing the Arrow topology coupling matrix.
    #         @param TZ = Transmission Zeros.
    #         @return MatQtriplet =  MatrixQ class instance, containing the Trisections topology coupling matrix. 
    #         
    #         [Cameron]
    #         
    def arrow2triplet(self, MatQarrow, TZ):
        Mtriplet = MatQarrow.copy()
        Mtriplet.setName('Triplet (N+2)')
        Mtriplet.setFileExt('.triplet')
        
        # Size of the matrix
        N = np.size(Mtriplet.M,0) - 2
        
        numTZ = len(TZ)
        for i in range(numTZ):
            N4 = i
            N3 = (N-3)-N4*2
            if N3>1:
                w01 = TZ[i]/1j
                rotAng = np.arctan(Mtriplet.M[N-1, N]/(w01+Mtriplet.M[N, N]))
                Mtriplet.rotateMatrix('trigonometric', N-1, N, rotAng, flagQ=False)
                for j in range(1, N3+1):
                    rotAng = Mtriplet.rotAngleEliminate('trigonometric', N-j-1, N-j, N-j-1, N-j+1)
                    Mtriplet.rotateMatrix('trigonometric', N-j-1, N-j, rotAng, flagQ=False)
            else:
                return None
            
        Mtriplet.Qresonators()
        return Mtriplet
        
        
    ##
    # 
    #         This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Cul de Sac coupling matrix \f$Mculdesac\f$.
    #         @param MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         @param TZ = Transmission Zeros.
    #         @return MatQculdesac = MatrixQ class instance, containing the Cul de Sac coupling matrix.
    #         
    #         [Cameron]
    #         
    def fcm2culdesac(self, MatQfcm, TZ):
    
        Mculdesac = MatQfcm.copy()
        Mculdesac.setName('Cul de Sac (N+2)')
        Mculdesac.setFileExt('.cds')
        
        # Size of the matrix
        N = np.size(Mculdesac.M,0) - 2
        
        numTZ = len(TZ)
        if N < 5 or numTZ > N-3:
            return None
        
        if N%2:
            R = (N-3)/2
            for r in range(1, R+1):
                i = (N+1)/2 - r
                j = (N+1)/2 + r
                rotAng = Mculdesac.rotAngleEliminate('trigonometric', i, j, i, j-1)
                Mculdesac.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        else:
            R = (N-2)/2
            for r in range(1, R+1):
                i = (N+2)/2 - r
                j = N/2 + r
                rotAng = Mculdesac.rotAngleEliminate('trigonometric', i, j, i, j)
                Mculdesac.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        
        Mculdesac.Qresonators()
        return Mculdesac
    
    ##
    # 
    #         This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Cascaded Quartets coupling matrix \f$Mcqs\f$.
    #         @param MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         @param r = If we want to shift the quartet, r indicates the position of the first resonator node of the quartet. It should be 0 if we do not want to shift it.
    #         @param UpDownNone = 'Up' : Shift the quartet above the diagonal. 'Down' : Shift the quartet below the diagonal. 'None' : No shifting.
    #         @return MatQcqs = MatrixQ class instance, containing the Cascaded Quartets coupling matrix.
    #         
    #         [Cameron]
    #         
    def fcm2cqs(self, MatQfcm, r, UpDownNone):
        
        Mcqs = MatQfcm.copy()
        Mcqs.setName('Cascaded Quartets (N+2)')
        Mcqs.setFileExt('.cqs')
        
        # Size of the matrix
        N = np.size(Mcqs.M,0) - 2
        
        if N < 8:
            return None
        
        R=(N-8)/2;

        # Four rotations are applied
        i = 3+R
        j = 5+R
        A = (Mcqs.M[2+R,7+R]*Mcqs.M[3+R,4+R]*Mcqs.M[4+R,5+R]-Mcqs.M[2+R,3+R]*Mcqs.M[5+R,6+R]*Mcqs.M[6+R,7+R]+Mcqs.M[2+R,7+R]*Mcqs.M[3+R,6+R]*Mcqs.M[5+R,6+R])
        B = (Mcqs.M[2+R,3+R]*Mcqs.M[3+R,6+R]*Mcqs.M[6+R,7+R]-Mcqs.M[2+R,7+R]*(Mcqs.M[3+R,4+R]**2-Mcqs.M[4+R,5+R]**2-Mcqs.M[5+R,6+R]**2+Mcqs.M[3+R,6+R]**2))
        C = -Mcqs.M[2+R,7+R]*(Mcqs.M[3+R,6+R]*Mcqs.M[5+R,6+R]+Mcqs.M[3+R,4+R]*Mcqs.M[4+R,5+R])
        rotAng = np.arctan((-B+np.sqrt(B**2-4*A*C))/(2*A))
        Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        
        i = 4+R
        j = 6+R
        rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, 3+R, j)
        Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        
        i = 5+R
        j = 7+R
        rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, 4+R, j)
        Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        
        i = 2+R
        j = 4+R
        rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, i, 5+R)
        Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        
        # If we want to shift the quartet (r is the value of the shift)
        if UpDownNone == 'Up':
            # up the diagonal
            i = r
            j = r+2
            rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, i, r+3)
            Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        elif UpDownNone == 'Down':
            # down the diagonal
            i = r+1
            j = r+3
            rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, r, j)
            Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        elif UpDownNone == 'None':
            None 
        else:
            raise synthError,  u"Incorrect command. UpDownNone must be 'Up', 'Down' or 'None'."
    
        Mcqs.Qresonators()
        return Mcqs
    
    ##
    # 
    #         This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Pfitzenmaier coupling matrix \f$Mpfi\f$.
    #         @param MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         @return MatQpfi = MatrixQ class instance, containing the Pfitzenmaier coupling matrix.
    #         
    #         [Cameron]
    #         
    def fcm2pfitzenmaier(self, MatQfcm):
    
        Mpfi = MatQfcm.copy()
        Mpfi.setName('Pfitzenmaier (N+2)')
        Mpfi.setFileExt('.pfi')
        
        # Size of the matrix
        N = np.size(Mpfi.M,0) - 2
        
        R = (N-4)/2
        for r in range(1, R+1):
            i = r+1
            j = N-i
            rotAng = Mpfi.rotAngleEliminate('trigonometric', i, j, i, N-r)
            Mpfi.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
    
        Mpfi.Qresonators()
        return Mpfi

    ##
    # 
    #         This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Inline Asymmetric coupling matrix \f$Minasym\f$.
    #         @param MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         @return MatQinasym = MatrixQ class instance, containing the Inline Asymmetric coupling matrix.
    #         
    #         [Cameron]
    #         
    def fcm2inlineAsymmetric(self, MatQfcm):
    
        Minasym = MatQfcm.copy()
        Minasym.setName('Inline Asymmetric (N+2)')
        Minasym.setFileExt('.ina')
        
        # Size of the matrix
        N = np.size(Minasym.M,0) - 2
        
        if N == 6:
            i = 2
            j = 4
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 5)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        elif N == 8:
            i = 4
            j = 6
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 2
            j = 4
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 7)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 3
            j = 5
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 2, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 5
            j = 7
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 4, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        elif N == 10:
            i = 4
            j = 6
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 7)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 6
            j = 8
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 7
            j = 9
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 6, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
        elif N == 12:
            i = 5
            j = 9
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 4, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 3
            j = 5
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 10)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 2
            j = 4
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 5)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 6
            j = 8
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 7
            j = 9
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 6, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 8
            j = 10
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 5, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
            i = 9
            j = 11
            rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 8, j)
            Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
    
        Minasym.Qresonators()
        return Minasym

    ##
    # 
    #         This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Inline Symmetric coupling matrix \f$Minsym\f$.
    #         @param MatQfcm =  MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix. 
    #         @param TZ = Transmission Zeros.
    #         @return MatQinsym = MatrixQ class instance, containing the Inline Symmetric coupling matrix.
    #         
    #         [Cameron]
    #         
    def fcm2inlineSymmetric(self, MatQfcm, TZ):
    
        Minsym = MatQfcm.copy()
        Minsym.setName('Inline Symmetric (N+2)')
        Minsym.setFileExt('.ins')
        
        # Size of the matrix
        N = np.size(Minsym.M,0) - 2
        
        numTZ = len(TZ)
        if (N == 6 and numTZ == 2) or ((N == 8 or N == 10) and numTZ == 4) or (N == 12 and numTZ == 6):
            None
        else:
            return None
        
        Meven1 = Minsym.copy()
        Meven1.M = Meven1.M[1:N/2+1, 1:N/2+1]
        Meven1.extraNodesL = Meven1.extraNodesL - 1
        Meven1.extraNodesS = Meven1.extraNodesS - 1
        Meven2 = Minsym.copy()
        Meven2.M = Meven2.M[1:N/2+1, N/2+1:N+1]
        Meven2.extraNodesL = Meven2.extraNodesL - 1
        Meven2.extraNodesS = Meven2.extraNodesS - 1
        Meven = Minsym.copy()
        Meven.M = Meven1.M + np.transpose(np.rot90(Meven2.M))
        Meven.extraNodesL = Meven.extraNodesL - 1
        Meven.extraNodesS = Meven.extraNodesS - 1
        
        if N == 6:
            i = 2
            j = 3
            rotAng = np.arctan((Minsym.M[2,3]-np.sqrt(((Minsym.M[2,3])**2)-(Minsym.M[2,5]*Minsym.M[3,4])))/Minsym.M[3,4])
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng, flagQ=False)
        elif N == 8:
            i = 3
            j = 4
            rotAng1 = np.arctan((Minsym.M[2,7]*Minsym.M[3,4]-np.sqrt((Minsym.M[2,7]**2)*(Minsym.M[3,4]**2)+(Minsym.M[2,7]*Minsym.M[4,5])*((Minsym.M[2,3]**2)-Minsym.M[2,7]*Minsym.M[3,6])))/(Minsym.M[2,3]**2-Minsym.M[2,7]*Minsym.M[3,6]))
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
            i = 2
            j = 4
            rotAng2 = np.arctan(Minsym.M[2,7]/(Minsym.M[2,3]*np.sin(rotAng1)))
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
        elif N == 10:
            i = 4
            j = 5
            rotAng1 = np.arctan((Minsym.M[4,5]-np.sqrt(Minsym.M[4,5]**2-Minsym.M[4,7]*Minsym.M[5,6]))/Minsym.M[5,6])
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
            i = 3
            j = 5
            rotAng2 = np.arctan((np.sin(rotAng1)*Minsym.M[3,4]-np.sqrt((np.sin(rotAng1)**2*Minsym.M[3,4]**2-Minsym.M[3,8]*(Minsym.M[4,7]+Minsym.M[5,6]))))/(Minsym.M[4,7]+Minsym.M[5,6]))
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
            i = 2
            j = 4
            rotAng3 = np.arctan(np.tan(rotAng2)*Minsym.M[2,3]/(Minsym.M[4,5]-np.tan(rotAng1)*Minsym.M[5,6]+np.cos(rotAng1)*np.tan(rotAng2)*Minsym.M[3,4]))
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng3, flagQ=False)
        elif N == 12:
            a0 = Minsym.M[5,8]
            a1 = 2*Minsym.M[4,5]
            a2 = Minsym.M[4,9]-((Minsym.M[3,4]**2)/Minsym.M[3,10])
            a3 = Minsym.M[6,7]
            a4 = -2*Minsym.M[5,6]
            b0 = Minsym.M[4,9]*Minsym.M[6,7]
            b1 = -2*Minsym.M[4,5]*Minsym.M[6,7]
            b2 = Minsym.M[5,8]*Minsym.M[6,7]-Minsym.M[5,6]**2
            b3 = -2*Minsym.M[5,6]*Minsym.M[4,5]
            b4 = Minsym.M[4,9]*Minsym.M[5,8]-Minsym.M[4,5]**2
            b5 = 2*Minsym.M[5,6]*Minsym.M[4,9]
            c=a3*b5-a4*b4
            c0 = (a0*b4-a3*b0)/c
            c1 = (a1*b4-a3*b1)/c
            c2 = (a2*b4-a3*b2)/c
            c3 = a3*b3/c
            d = a2*c3**2+a3*c2**2
            d0 = (a0+a3*c0**2+a4*c0)/d
            d1 = (a1+2*a0*c3+2*a3*c0*c1+a4*(c1+c0*c3))/d 
            d2 = (a2+2*a1*c3+a0*c3**2+a3*(c1**2+2*c0*c2)+a4*(c2+c1*c3))/d
            d3 = (2*a2*c3+a1*c3**2+2*a3*c1*c2+a4*c2*c3)/d
            t1 = np.roots([1., d3, d2, d1, d0])[0]
            i = 4
            j = 5
            rotAng1 = np.arctan(t1)
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
            i = 5
            j = 6
            rotAng2 = np.arctan((c0+c1*t1+c2*t1**2)/((1+c3*t1)*np.sqrt(1+t1**2)))
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
            i = 4
            j = 6
            rotAng3 = np.arctan(Meven.M[i-1, j-1]/Meven.M[j-1, j-1])
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng3, flagQ=False)
            i = 3
            j = 5
            rotAng4 = np.arctan(Meven.M[i-1, j-1]/Meven.M[j-1, j-1])
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng4, flagQ=False)
            i = 2
            j = 4
            rotAng5 = np.arctan(Meven.M[i-1, j]/Meven.M[j-1, j])
            Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng5, flagQ=False)
        
        Minline = np.zeros([N/2,N], dtype='complex128')
        # Upper diagonal
        for i in range(N/2):
            if i<N/2-1:
                Minline[i,i+1] = Meven.M[i,i+1]
                Minline[i+1,i] = Minline[i,i+1]
            else:
                Minline[i,i+1] = Meven.M[i,i]       

        # Other elements
        if N == 6:
            Minline[0,3] = Meven.M[0,2]
            Minline[2,5] = Minline[0,3]
        elif N == 8:
            Minline[0,3] = Meven.M[0,3]
            Minline[3,0] = Minline[0,3]
            Minline[2,5] = Meven.M[2,2]
        elif N == 10:
            Minline[0,3] = Meven.M[0,3]
            Minline[3,0] = Minline[0,3]
            Minline[2,5] = Meven.M[2,4]
            Minline[4,7] = Minline[2,5]
        elif N == 12:
            Minline[0,3] = Meven.M[0,3]
            Minline[3,0] = Minline[0,3]
            Minline[2,5] = Meven.M[2,5]
            Minline[4,7] = Meven.M[4,4]
            Minline[5,2] = Minline[2,5]

        # Fulfill the two lower quadrants of the matrix, which are symmetric to the upper ones
        Minline1=np.rot90(Minline)
        Minline2=np.rot90(Minline1)
        Minsym.M[1:N+1, 1:N+1]=np.concatenate([Minline, Minline2])

        Minsym.Qresonators()
        return Minsym
    
    ##
    # 
    #         Read coupling matrix from disk file.
    #         @param fileName = File name (string).
    #         @return MatQ = Coupling matrix (MatrixQ class instance.)
    #         
    def readCouplingMatrix(self, fileName):
        try:
            file = codecs.open(fileName,  'r',  'utf-8')
        except IOError,  err:
            raise parseError,  'I/O error: %s' % err
        stList = file.readlines()
        file.close()
        
        validLines = [ st for st in stList if not st.startswith('#') and not st.isspace() ] 

        pat = re.compile('^Name:(.+)')
        name = pat.match(validLines[0]).groups(1)[0].lstrip().rstrip()

        pat = re.compile('^extraNodesS:(.+)')
        extraNodesS = int(pat.match(validLines[1]).groups(1)[0])

        pat = re.compile('^extraNodesL:(.+)')
        extraNodesL = int(pat.match(validLines[2]).groups(1)[0])

        pat = re.compile('^nodeScalingFactor:(.+)')
        nodeScalingFactor = np.array( eval(pat.match(validLines[3]).groups(1)[0].strip()) )

        pat = re.compile('^flagRotateAllowed:(.+)')
        flagRotateAllowed = eval(pat.match(validLines[4]).groups(1)[0].strip())
        
        # Matrix order
        N = len(nodeScalingFactor)

        Mlist = []
        for line in validLines[6:6+N]:
            Mlist.append([complex(elem) for elem in line.split()])
        M = np.array(Mlist)
        
        MatQ = MatrixQ(self, M,  name,  extraNodesS,  extraNodesL,  self.FT,  nodeScalingFactor=nodeScalingFactor,  flagRotateAllowed=flagRotateAllowed)
        myPrint("Coupling matrix '%s' read from file: %s" % (name, os.path.basename(fileName)) )
        return MatQ
        

    ##
    # 
    #         Save coupling matrices in output file.
    #         The matrices are attributes of CouplingMatrices class, and of type numpy.array.
    #         @param P = Parameter class instance.
    #         @param SP = SParameters class instance.
    #         @param precCM   = Number of significant digits to save in coupling matrix and Q (int).
    #         @param precEP   = Number of significant digits to save in energy and power (int).
    #         
    def saveCouplingMatrices(self, P, SP, precCM, precEP):
        # Transversal coupling matrix
        myPrint('Coupling matrices:')
        for MatQ in self.listM:
            if MatQ.fileExt is None: continue
            myPrint(MatQ.name + ':')
            st = "   Q of resonators = %s" %  ["%.4g" % elem if abs(elem) != np.inf else "%s" % str(elem) for elem in MatQ.Q] # ["%.4g" % elem for elem in MatQ.Q]
            myPrint(st.replace("'", ""))
            fileName = P.outDirName + MatQ.fileExt
            MatQ.saveMat(fileName,  P, SP, precCM, precEP)
            myPrint("   Matrix saved to file: %s" % fileName)
 

    ##
    # 
    #         This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 4 into an N+4 matrix so that all the resonant nodes have equal quality factor Q. 
    #         (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
    #         @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
    #         
    #         For a better description see TN102.3 sec. 3.2.
    #         
    #         First of all we add two non-resonant nodes, what means converting the N+2 folded matrix into the N+4. Then we apply two hyperbolic rotations so that we have a uniform Q. For doing so, it can be proved that the necessary angle is (TN102.3 eq. (10)):
    #         \f[
    #         \alpha = \frac{1}{2}ln\left[\frac{\frac{M_{01}^{2}}{G}\pm \sqrt{4G_{1}\frac{M_{01}^{2}}{G}-4G_{1}^{2}+16M_{12}^{2}}}{2G_{1}+\frac{M_{01}^{2}}{G}+4M_{12}}\right]
    #         \f]
    #         Finally, we need to re-scale with a value \f$h\f$ so that the non-resonant nodes have a sum of imaginary parts equal to 0.
    #         
    def uniformQFCMOrder4(self, MatQfcm2):
        MatQ = MatQfcm2.copy()
        MatQ.setName('Folded with uniform Q (N+4)')
        MatQ.setFileExt('.fcm_uniQ')
        
        # Adition of the non-resonant nodes (N+4 matrix)
        MatQ.inflateMatrix('both')
        
        # A = -G1 = Im(Mfcm4_22)
        A=np.imag(MatQ.M[2, 2])
        # B = 2M_12= 2Mfcm4_23
        B=2*MatQ.M[2, 3]
        # C = - (M_01)^2 / G = (Mfcm4_12)^2 / Im(Mfcm4_11)
        C=MatQ.M[1, 2]**2 /np.imag(MatQ.M[1, 1])
        
        # rotation angle alpha from TN.102.3 eq. (10)
        alpha=0.5*np.log(1./2/(-2*A-C+2*B)*(-2.*C+4.*np.sqrt(-A**2-C*A+B**2)))
        
        # Aplication of the rotations to make Q uniform
        i=2;j=3;
        MatQ.rotateMatrix('hyperbolic',i,j,alpha, flagQ=False)
        i=5;j=4;
        MatQ.rotateMatrix('hyperbolic',i,j,alpha, flagQ=False)
        
        # Node Scaling to make the sum of the imaginary parts of the non-resonant nodes equal to zero, and therefore Q=0 in those nodes
        h=-np.imag(MatQ.M[3, 1])/np.imag(MatQ.M[1, 1])
        MatQ.scaleNode(1, h, flagQ=False)
        h=-np.imag(MatQ.M[4, 6])/np.imag(MatQ.M[6, 6])
        MatQ.scaleNode(6, h, flagQ=True)
        
        return MatQ

    ##
    # 
    #         This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 6 into an N+4 matrix so that all the resonant nodes have equal quality factor Q. 
    #         (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
    #         @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
    #         
    #         See TN102.3 sec. 3.3.        
    #         
    def uniformQFCMOrder6(self, MatQfcm2):
        MatQ = MatQfcm2.copy()
        MatQ.setName('Folded with uniform Q (N+4)')
        MatQ.setFileExt('.fcm_uniQ')
        
        # Obtain the rotation angle which makes Q uniform, with an optimization routine
        print "\nUNIFORM Q OPTIMIZATION:"
        angRotOpt = scop.fmin(self.uniformQOrder6Optimize, np.array([0,0]), args=(MatQ, 0), xtol=1.e-8, ftol=1.e-8, maxiter=20000, maxfun=20000, full_output=0, disp=True, retall=0, callback=None)
        print ""
        
        # Resulting matrix for the optimal angle found
        MatQ = self.uniformQOrder6Optimize(angRotOpt, MatQ, 1)
        
        # Conversion to the N+4 matrix taking into account the node scaling done k
        MatQ.inflateMatrix('both')
        
        return MatQ
        
    
    ##
    # 
    #         This function performs the hyperbolic rotations in TN102.3 sec.3.3.
    #         Given two angles angRot it performs the hyperbolic rotations with those angles 
    #         and sees the differences between the quality factors of the different nodes [TN102.3 sec.3.3].
    #         @param angRot = Vector containing two angles. The first angle corresponds with \f$\alpha\f$ in TN102.3 eq. (14) and the second one with \f$\alpha_{1}\f$
    #         @param MatQfcm2a = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         @param flag = If False, return goal; if True, return MatQfcm2.
    #         @return goal = Function we are minimizing. It contains the summation of the absolute values of the diferences in losses in the different nodes.
    #         @return MatQfcm2 = MatrixQ class instance, containing the transformed Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         
    def uniformQOrder6Optimize(self, angRot, MatQfcm2a, flag):
        MatQfcm2 = MatQfcm2a.copy()
        angle12 = angle65 = angRot[0]
        angle23 = angle54 = angRot[1]
        
        # Application of the rotations
        i=1; j=2; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle12, flagQ=False)
        i=6; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle65, flagQ=False)
        i=2; j=3; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle23, flagQ=False)
        i=5; j=4; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle54, flagQ=False)
    
        #Node scaling
        k=-np.imag(MatQfcm2.M[2,0])/np.imag(MatQfcm2.M[0,0])
        MatQfcm2.scaleNode( 0, k, flagQ=False)
        MatQfcm2.scaleNode( 7, k, flagQ=True)
        
        Q = np.array(MatQfcm2.Q)
        # Function to minimize
        goal = abs(1/Q[0]-1/Q[1]) + abs(1/Q[1]-1/Q[2]) + abs(1/Q[2]-1/Q[3]) + abs(1/Q[3]-1/Q[4]) + abs(1/Q[4]-1/Q[5])
        
        if flag:
            return MatQfcm2
        else:
            return goal
            
            
    ##
    # 
    #         This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 6 into an N+4 matrix so that all the resonant nodes have equal quality factor Q. 
    #         (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
    #         @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
    #         
    def uniformQFCMOrder6Case3(self, MatQfcm2):
        MatQ = MatQfcm2.copy()
        MatQ.setName('Folded with uniform Q (N+4)')
        MatQ.setFileExt('.fcm_uniQ')
        
        # Obtain the rotation angle which makes Q uniform, with an optimization routine
        print "\nUNIFORM Q OPTIMIZATION:"
        angRotOpt = scop.fmin_slsqp(self.uniformQOrder6Case3Optimize, np.array([0,0, 0, 0]) , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None,bounds = [], 
                                                        fprime = None, fprime_eqcons=None, fprime_ieqcons=None, args = (MatQ, False), iter = 20000, acc = 1.0E-10, iprint = 1, full_output = 0)
        print ""
        
        # Resulting matrix for the optimal angle found
        MatQ = self.uniformQOrder6Case3Optimize(angRotOpt, MatQ, True)
        
        # Conversion to the N+4 matrix taking into account the node scaling done k
        MatQ.inflateMatrix('both')
        return MatQ
    

    ##
    # 
    #         This function performs the hyperbolic rotations in TN102.3.
    #         Given four angles angRot it performs the hyperbolic rotations with those angles.
    #         and sees the differences between the quality factors of the different nodes [TN102.3].
    #         @param angRot = Vector containing two angles. The first angle corresponds with \f$\alpha\f$ in TN102.3 eq. (14) and the second one with \f$\alpha_{1}\f$
    #         @param MatQfcm2a = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         @param flag = If False, return goal; if True, return MatQfcm2.
    #         @return goal = Function we are minimizing. It contains the summation of the absolute values of the diferences in losses in the different nodes.
    #         @return MatQfcm2 = MatrixQ class instance, containing the transformed Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix. 
    #         
    def uniformQOrder6Case3Optimize(self, angRot, MatQfcm2a, flag):
        MatQfcm2 = MatQfcm2a.copy()
        
        angle14 = angRot[0];
        angle23 = angRot[1];
        angle12 = angRot[2];
        angle34 = angRot[3];
        
        # Application of the rotations
        i=2; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle14, flagQ=False)
        i=3; j=4; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle23, flagQ=False)
        i=2; j=3; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle12, flagQ=False)
        i=4; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle34, flagQ=True)
        
        Q = MatQfcm2.Q
        # Function to minimize: Q of the four middle resonators, not the 1st and the last ones.
        goal = abs(1/Q[1]-1/Q[2]) + abs(1/Q[2]-1/Q[3]) + abs(1/Q[3]-1/Q[4])
        
        if flag:
            return  MatQfcm2
        else:
            return goal

# End class: CouplingMatrices
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin: Frequency Transformation classes

##
# 
#     Frequency transformation fro band pass to low pass prototype.
#     
class FrequencyTransformBP(object):
    ##
    # 
    #         Create a frequency transformation instance.
    #         
    def __init__(self, P):
        self.BW = P.BW
        self.f0 = P.f0
        
        ##
        # Fractional bandwidth.
        self.FBW = P.BW/P.f0
        
        
    ##
    # 
    #         Normalize frequency to lowpass prototype.
    #         @param freq = Vector with frequency values to convert (in Hz). Numpy array.
    #         @return Vector with lowpass prototype normalized frequency values. Numpy array.
    # 
    #         [Cameron eq. (3.113)]:
    #         
    #         \f$ \omega ' = \frac{f_0}{\Delta f} \left( \frac{f}{f_0} - \frac{f_0}{f} \right) \f$
    #         
    def normFreq(self,  freq):
#            # Obsolete code. It can be removed unless proved otherwise
#        if type(freq) != types.ListType: freq = [ freq ]
#        ret_freq = [ (self.f0 / self.BW) * (f/self.f0 - self.f0/f) if f is not None else None for f in freq ]
#        if len(ret_freq) == 1: ret_freq = ret_freq[0]
        
        return (self.f0 / self.BW) * (freq/self.f0 - self.f0/freq)

    
    ##
    # 
    #         Unnormalize frequency from lowpass prototype.
    #         @param w = Vector with normalized frequency values to convert. Numpy array.
    #         @return Vector with unnormalized frequency values (in Hz). Numpy array.
    # 
    #         \f[
    #             \omega ' = \frac{f_0}{\Delta f} \left( \frac{f}{f_0} - \frac{f_0}{f} \right) 
    #         \f]
    #         
    #         \f[
    #         f^2 - \omega ' \Delta f f - f_0^2 = 0 
    #         \f]
    #         
    #         \f[
    #         f = \frac {\omega '  \Delta f + \sqrt{( \omega '  \Delta f)^2 + 4 f_0^2 } } {2} 
    #         \f]
    #         with \f$ \omega'' = \omega' \Delta f / f_0 \f$ :
    #         \f[
    #               f = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
    #         \f]
    #         
    def unormFreq(self,  w):
#        b = -w*self.BW
#        return (-b+np.sqrt(b**2+4*self.f0**2))/2
        wp = w*self.FBW
        return (self.f0/2) * (wp + np.sqrt(wp**2 +4)) 
    
        
    ##
    # 
    #         Convert lowpass prototype group delay to unnormalized group delay.
    #         @param freq = Vector with unnormalized frequency values (in Hz). Numpy array.
    #         @param groupDelay = Vector with lowpass prototype group delay values to convert. Numpy array.
    #         @return Vector with unnormalized filter group delay values (in seconds). Numpy array.
    #         
    #         [Cameron, eq. (3.90)]:
    #         
    #         \f$ \tau_{BPF}= \frac{\omega_0}{\Delta \omega} \left( \frac{1}{\omega_0} -  \frac{\omega_0}{\omega^2}  \right) \tau ' = \frac{1}{2 \pi \Delta f}  \left( 1 +  \left( \frac{f_0}{f} \right)^2  \right) \tau '  \f$
    #         
    def unormGroupDelay(self,  freq,  groupDelay):
        return groupDelay * (1 + (self.f0/freq)**2) / (2*pi*self.BW) 
        
    ##
    # 
    #         Convert unnormalized group delay to lowpass prototype group delay.
    #         @param groupDelay = Double with unnormalized group delay (in seconds) value to convert. 
    #         @return groupDelay = normalized filter group delay value.
    #         
    #         [Cameron, eq. (3.91)]:
    #         
    #         \f$ \tau_{0}= \frac{2}{\Delta \omega}  \tau_{0} '  \f$
    #         
    def normGroupDelay(self,  groupDelay):
        return groupDelay * (2*pi*self.BW) / 2 
        
# End:  Frequency Transformation classes
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#